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Abstract 

The  non-linear  dynamic  behavior  of  suddenly  loaded  composite  shells  is  the  subject  of  this 
research.  The  objective  of  the  research  is  two-fold:  (1)  to  characterize  the  apparently  random 
physical  behavior  sometimes  observed  in  the  finite  element  analysis  as  either  numerical  instability, 
chaotic  behavior,  or  both;  and  (2)  to  develop  criteria  for  choosing  time  steps  for  the  finite  element 
model,  referred  to  as  DSHELL.  This  displacement  based  finite  element  code  is  capable  of  dynamic 
analysis  of  plates,  arches,  and  cylindrical  shells  undergoing  large  displacements  and  moderately 
large  rotations  during  deformation.  DSHELL  uses  a  36-DOF  isoparametric  shell  element  to  obtain 
numerical  results.  A  simplified  shell  model  (MSHELL)  is  developed  to  provide  a  shell  analog  that 
requires  much  less  computer  time  to  run  than  the  full-up  finite  element  model.  The  results  of  the 
investigations  using  this  simplified  model  are  then  applied  toward  understanding  behavior  seen  in 
the  finite  element  code.  The  simplified  model  and  the  composite  shell  exhibit  chaotic  behavior  after 
collapse.  For  loads  not  sufficient  to  cause  collapse,  the  composite  shell  experiences  oscillations  at 
the  panel  edges  of  a  type  associated  with  near-chaotic  behavior.  Significant  in-plane  oscillations 
(parametric  resonances)  are  exhibited  by  the  suddenly-loaded  composite  shell.  Extensive  use  of 
the  Fast  Fourier  Transform  (FFT)  is  made  to  investigate  behavior  of  the  MSHELL  and  DSHELL 
models  in  the  frequency  domain.  The  composite  shell  is  compared  to  an  isotropic  shell  of  identical 
physical  dimensions  and  static  collapse  strength,  showing  the  primary  effect  of  isotropy  to  be 
increased  flexibility  in  all  displacement  directions.  A  method  is  developed  using  MSHELL  for 
choosing  appropriate  time  steps  for  analyses  using  DSHELL. 
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The  Non-Linear  Dynamics  of  Composite  Cylindrical  Shells 


I.  INTRODUCTION 

The  use  of  laminated  composite  materials  in  structures  continues  to  grow.  The  high  strength- 
to-weight  ratios  offered  by  these  materials  makes  them  ideally  suited  to  cierospace  applications. 
One  of  the  ways  laminated  fiber  composites  achieve  this  property  is  to  orient  the  fibers  in  preferred 
directions  depending  upon  the  required  strength  of  the  material  and  its  load-bearing  orientation. 
This  technique,  however,  often  renders  the  optimized  composite  substantially  weaker  in  the  non¬ 
preferred  directions  than  its  metallic  counterpart,  making  dynamic  stability  an  issue  (9:1). 

Most  research  involving  the  study  of  shell  structures  from  a  dynamic  point  of  view  exploits  the 
shell’s  definition — it  is  thin.  That  is,  the  thickness  dimension  is  much  smaller  than  the  shell’s  other 
dimensions.  In  many  cases,  this  assumption  allows  the  analyst  to  treat  some  quantities  as  negligible 
compared  to  others.  In  these  approaches,  the  in-plane  activity  of  the  shell  (displacements,  strains, 
and  stresses)  are  assumed  to  be  much  more  important  than  through-the-thickness  phenomena. 
The  outgrowth  of  this  approach  is  a  method  of  describing  the  behavior  of  the  shell  by  reference  to 
the  behavior  of  its  midsurface  (or  datum  surface).  In  many  cases,  particularly  those  dealing  with 
small  displacements  and/or  isotropic  materials,  this  method  is  sufficient.  However,  the  analysis  of 
laminated  composite  materials,  where  material  properties  may  vary  drastically  through  the  shell  s 
thickness  dimension,  requires  careful  evaluation  of  which  quantities  may  or  may  not  be  neglected . 

1 . 1  Nonlinear  Dynamics 

If  the  analyst  is  interested  in  postbuckling  behavior,  simple  linear  bifurcation  buckling  analy¬ 
ses  will  not  suffice,  and  solution  of  the  nonlinear  equations  resulting  from  a  collapse  analysis  will  be 
required  (9:2).  Axially  loaded  columns  and  edge  loaded  flat  plates  both  exhibit  instabilities  char¬ 
acterized  by  very  small  deformation  prior  to  buckling.  On  the  other  hand,  transversely  loaded  shell 
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structures  may  change  shape  and  stress  distribution  significantly  prior  to  collapse.  Attempting  to 
find  the  critical  collapse  load  of  shells  through  linear  bifurcation  techniques  are  often  significantly 
in  error  (19:234). 

Through  linear  eigenvalue  analysis  using  equilibrium  or  energy  methods  the  critical  collapse 
load,  Pc  ,  of  many  simple  structures  may  be  solved  for  in  closed  form.  However,  the  calculated 
value  of  Pc  may  be  very  much  at  odds  with  experimental  values  which  are  often  only  1/3  as  much 
as  the  calculated  value.  Furthermore,  the  experimental  data  frequently  have  large  amounts  of 
scatter.  For  some  30  years,  the  discrepancy  between  experimental  and  calculated  values  of  Pc 
baffled  investigators.  The  source  of  the  discrepancy  was  found  to  be  the  nature  of  the  initial 
postbuckling  branching.  If  the  postbuckling  branching  was  stable  (i.e.,  increasing  the  displacement 
required  increasing  the  load)  then  the  calculated  and  experimental  data  agreed  fairly  well.  If  the 
postbuckling  branching  was  unstable  (i.e.,  increasing  the  displacement  required  decreasing  the  load) 
then  the  calculated  and  experimental  data  were  frequently  at  odds.  (10:409-413) 

Figure  3.1  on  page  3-2  illustrates  the  type  of  behavior  seen  in  unstable  buckling.  In  those 
portions  of  the  curve  where  the  slope  dP/dw  is  negative,  a  decrease  in  load,  P  ,  is  required  to  reach 
an  equilibrium  state  at  increased  deflection.  It  is  for  this  reason  that  buckling  shells  and  shell-like 
structures  (like  the  simplified  model  used  in  the  present  research)  are  not  correctly  analyzed  using 
linear  bifurcation  techniques. 

1.2  Chaos 

Another  interesting  feature  of  the  dynamic  analysis  of  shells  has  been  the  appearance  of 
chaoslike  phenomena.  Chaos  is  the  term  used  to  define  the  behavior  of  dynamical  systems  whose 
motion  cannot  be  predicted  far  into  the  future  (17:1).  There  is  a  difference,  however,  between 
so-called  random  and  chaotic  motions.  The  former  term  is  reserved  for  cases  where  parameters 
may  be  ill-defined  or  unknown  except  in  perhaps  a  statistical  sense.  The  term  chaotic  is  reserved 
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for  those  cases  where  the  there  are  no  unpredictable  inputs  or  system  parameters  (17:3).  Chaotic 
motion  includes  the  following  chuacteristics  (17:262): 

1.  The  motion  is  sensitive  to  changes  in  initial  conditions 

2.  Trajectories  starting  from  slightly  different  initial  conditions  diverge  exponentially 

3.  The  Lyapunov  exponent  for  the  motion  is  positive 

Many  physical  systems  exhibit  chaotic  vibrations,  including  the  vibrations  of  buckled  elastic  struc¬ 
tures  (17:9). 

The  nonlinear  equations  of  equilibrium  derived  for  the  DSHELL  finite  element  model  have 
both  of  the  necessary  ingredients  for  chaos:  (a)  at  least  three  independent  dynamical  variables 
(DSHELL  has  seven)  and  (b)  nonlinear  terms  in  the  coupled  equations  (there  are  many).  Meeting 
these  criteria  admits  the  following  behavior  (2:3); 

•  trajectories  in  phase  space  may  diverge 

•  confinement  (boundedness)  of  motion  to  a  finite  region  of  the  phase  space 

•  uniqueness  of  the  trajectory,  with  the  added  feature  that  two  trajectories  corresponding  to 
similar  energies  will  pass  very  close  to  each  other,  but  the  orbits  will  not  cross  each  other. 

Taylor  (28)  and  others  have  witnessed  an  apparently  random  behavior  in  both  pre-  and  post-buckled 
states,  and  the  present  research  shows  that  there  are  two  sources  of  this  behavior  —  one  of  them 
being  chaoslike.  The  deformed  geometries  of  the  shell  structures  differ  markedly  depending  upon 
whether  the  random  behavior  is  due  to  chaos  or  numerical  instability. 
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l.S  Analytical  Method 


An  analysis  of  the  laminated  composite  shell  subjected  to  a  sudden  transverse  load  could  take 
many  forms  depending  on  how  the  problem  is  cast.  Table  1.1  outlines  some  of  the  problem  features 
and  the  associated  considerations. 


Feature 

Consideration 

magnitude  of  displacements  and  rotations 

geometric  linearity 

magnitude  of  strains;  material  properties 

material  linearity;  plasticity 

rates  of  strain;  time  dependency 

creep  and  stress  rela.xation  (viscoelasticity) 

temperature;  temperature  changes 

thermal  effects 

Table  1.1.  Features  and  considerations  for  analysis 


i 

I 

Of  the  considerations  in  Table  1.1,  only  geometric  nonlinearity  is  addressed  in  the  present 
research.  Displacements  and  rotations  are  large,  so  even  if  small  strains  allow  the  assumption  of 
material  linearity,  geometric  linearity  cannot  be  claimed. 

The  variation  of  the  total  potential  energy  is  used  to  develop  the  governing  equations  for 
the  laminated  composite  shell.  The  result  of  applying  this  technique  is  five  coupled  nonlinear 
partial  differential  equations  that  describe  the  shell’s  equilibrium.  The  finite  element  method  is  a 
powerful  numeric^d  tool  for  solving  such  equations,  and  hence  is  the  method  of  choice  for  the  current 
research  (9:78).  The  finite  element  code  is  computationally  intensive,  often  modelling  hundreds  of 
nonlinearly  interacting  degrees  of  freedom  (DOF).  The  matrix  manipulations  require  large  amounts 
of  computational  storage  and  time.  For  this  reason,  the  MSHELL  .Nimplified  model  is  developed 
as  an  analog  to  the  finite  element  model.  The  MSHELL  simplified  model  has  many  of  the  shell’s 
features,  but  requires  comparatively  little  computational  time.  Hence,  long  simulations  may  be 
performed  using  MSHELL  in  an  attempt  to  characterize  the  phenomena  seen  in  the  finite  element 
code. 
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1-4  Previous  Work 


The  development  of  shell  theory  has  its  roots  in  plate  theory,  the  plate  being  a  special  case  of 
the  shell  (a  shell  with  infinite  radius  of  curvature  in  all  in-plane  directions).  Many  theories  have  been 
developed  to  treat  the  behavior  of  plates  and  shells.  Each  has  certain  advantages  and  drawbacks, 
and  each  has  some  limited  range  of  applicability.  The  most  general  equations  of  elasticity  governing 
the  behavior  of  plates  and  shells  are  cumbersome,  and  for  this  recison  each  theory  trades  at  me 
point  on  applicability  (generality)  versus  ease  of  use. 

Recent  work  in  the  area  of  shell  dynamics  is  wide-ranging.  Siniitses  (26)  offers  an  elegant 
analytical  treatment  of  thin  laminated  shells  (wherein  each  lamina  is  orthotropici  subjected  to 
sudden  loads.  Based  upon  his  general  assumptions,  In-  derives  the  field  equations  bcised  upon 
Donnell-type  kinematic  relations.  He  assumes  normals  to  the  shell  surface  remain  normal  after 
deformation,  thus  neglecting  through  the  thickness  shear,  and  further  restricts  the  problem  to  small 
strains  and  small  rotations.  He  presents  a  closed  form  static  solution  methodology  as  well  as  a  finite 
difference  dynamic  solution  methodology  for  the  thin  shell.  Simitses  also  describes  the  phenomenon 
of  parametric  resonance,  in  which  the  forcing  function  causes  “. . .  a  resonance  in  a  secondary  mode 
of  vibration  or  a  mode  other  than  the  one  directly  excited  by  its  presence.”  (26:260).  This  is  seen 
in  the  transversely  loaded  cylindrical  shell,  which  not  only  oscillates  in  the  tr8uisverse  direction, 
but  experiences  in-plane  oscillations  as  well. 

A  perturbation  procedure  is  used  by  Raouf  and  Palazotto  (21)  to  derive  a  set  of  asymptot¬ 
ically  consistent  nonlinear  equations  of  motion  describing  the  nonlinear  dynamics  of  an  arbitrar¬ 
ily  laminated  composite  cylindrical  shell  in  cylindrical  bending.  Unlike  the  present  research,  in 
which  transverse  shear  strains  are  allowed  while  strain  in  the  thickness  direction  is  not,  Raouf  and 
Palazotto  prohibit  transverse  shear  strains  but  allow  for  slight  compressibility  of  the  panel.  The 
three-dimensional  shell  equations  of  elasticity  are  reduced  to  two-dimensional  field  equations  for 
the  panel  using  an  order  of  magnitude  study.  The  analysis  reveals  a  threshold  value  for  the  ampli- 


tude  of  excitation  that  causes  saturation  of  the  amplitude  of  the  directly  excited  mode,  resulting 
in  a  nonlinear  response  of  the  coupled  mode  which  eventually  dominates  the  response.  In  subse¬ 
quent  work,  Raouf  and  Palazotto  (20)  take  the  quadratic  nonlinearities  used  in  the  perturbation 
procedure  and  perform  a  single-mode  spatial  discretization  to  derive  a  nonlinear  equation  with 
quadratic  and  cubic  nonlinearities.  The  results  of  this  study  of  a  symmetrically  laminated  eight-ply 
graphite/epoxy  cylindrical  shell  panel  show  that  nonlinearities  of  the  hardening  type  are  present 
even  for  small-amplitude  oscillations.  It  had  been  thought  that  linear  analysis  was  sufficient  for 
such  cases. 

Palazotto  and  Dennis  (19)  and  Reddy  and  Liu  (22)  have  both  presented  higher  order  theories 
incorporating  transverse  shear  effects,  and  this  is  certainly  a  key  feature  of  the  DSHELL  model  used 
in  the  current  research.  Incorporating  transverse  shear  effects  removes  the  limits  on  the  magnitudes 
of  deformations  required  by  first-order  theory.  Tsai  and  Palazotto  (30)  compared  the  results  of  this 
36  degree-of-freedom  thin  shell  element  based  code  with  analyses  from  other  sources,  demonstrating 
the  validity  and  accuracy  of  the  model. 

Smith  (27)  investigated  eight  variations  of  higher-order  theory  for  cylindrical  shells.  He 
showed  that  while  higher-order  effects  are  significant  during  the  collapse  phase  of  the  static  load- 
displacement  equilibrium,  the  extra  computational  burden  imposed  is  excessive  for  typical  engineer¬ 
ing  problems.  The  SHELL  code  uses  only  the  linear  terms  of  the  transverse  shear  strain  equations, 
but  uses  cubic  displacement  equations  (see  Section  2.1.1). 

Silva  used  SHELL  to  investigate  (statically)  the  behavior  of  cross-  and  angle-ply  laminates 
with  different  boundary  conditions,  evaluating  some  48  different  cases.  He  demonstrates  that  trans¬ 
verse  shear  has  a  significant  effect  on  the  deformation  of  composite  shells  subjected  to  a  transverse 
load,  and  that  the  angle-ply  laminate  developed  more  transverse  shear  strain  than  the  cross-ply. 
Silva  also  used  the  Donnell  shell  equations  to  show  how  they  develop  significant  inaccuracy  when 
the  rotations  of  the  midsurface  exceed  about  15®.  He  also  shows  that  the  angle-ply  laminate  can- 
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not  be  successfully  modeled  using  symmetry  because  the  twisting  generated  by  the  45“  plies  is  not 
captured  (25). 

Palazotto,  et  al.  (18)  used  DSHELL  to  study  dynamic  behavior  of  composite  arches  and 
shells,  demonstrating  the  importance  of  dynamic  analysis  (rather  than  static)  in  exploring  the  shell 
buckling  phenomenon.  Shell  depth  and  effects  of  anisotropy  were  examined  in  analyzing  the  shell 
snapping  phenomenon,  during  which  energy  is  transferred  from  the  potential  to  the  kinetic  as  the 
applied  load  reaches  a  critical  value. 

Chien  and  Palazotto  (6)  investigated  dynamic  buckling  of  composite  cylindrical  shells  sub¬ 
jected  to  a  step  loading.  They  found  that  the  dynamic  buckling  load  of  the  model  could  be  found 
through  an  evaluation  of  the  motion  at  the  node  under  the  transverse  load.  For  a  multi-degree-of- 
freedom  system,  the  existence  of  at  least  one  eigenvalue  with  a  positive  real  part  in  the  Jacobian 
matrix  associated  with  the  Hamiltonian  equation  indicates  dynamic  instability,  but  Chien  and 
Palazotto  found  that  motion  at  a  single  degree  of  freedom  could  be  evaluated  for  several  step 
loads  to  develop  a  “dynamic  equilibrium  curve”  that  could  predict  the  dynamic  snap  (or  collapse) 
load  (6:723).  This  procedure  is  used  cautiously,  and  other  arbitrary  degrees  of  freedom  are  checked 
for  the  potential  occurrence  of  parametric  resonance.  The  equation  defining  dynamic  instability  is 
shown  to  be 

Wp  =  Wp  =  0  (1.1) 

where  tip  and  tip  represent  the  velocity  and  acceleration  at  the  degree  of  freedom  in  question.  They 
show  that  the  presence  of  an  inflection  point  in  the  displacement  versus  time  curve  is  equivalent  to 
this  dynamic  buckling  criterion.  Cross-ply  laminates  are  investigated  in  the  form  of  an  arch  and  a 
cylindrical  shell,  eis  is  the  isotropic  panel. 

Maestrello,  et  al.  ( 16),  found  panels  of  aluminum  and  composite  material  vibrated  chaotically 
when  subjected  to  acoustic  excitation.  By  varying  the  sound  pressure  level,  they  were  able  to 
observe  the  transition  from  linear  to  nonlinear  to  chaotic  vibrations.  In  the  tran.sition  from  linear 
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to  nonlinear  vibration  they  note  the  appearance  of  many  sub-harmonics  in  the  frequency  response  of 
the  panel.  A  method  formulated  by  Zeng,  et  al.  (32)  was  used  to  estimate  the  Lyapunov  exponents 
from  their  limited  data,  and  the  results  indicated  very  chaotic  behavior  under  certain  conditions. 

Blair,  et  al.  (4)  used  a  two  rigid-link  model  to  explore  the  nonlinear  dynamic  response  of 
shallow  arches,  finding  regions  of  chaotic  motion.  They  used  the  method  of  harmonic  balance, 
coupled  with  a  continuation  scheme  to  find  solutions  to  a  wide  range  of  externally  applied  forcing 
functions. 

Currently  at  AFIT,  Captain  Scott  Schimmels  is  developing  codes  to  incorporate  plasticity  in 
the  finite  element  analysis  of  plates,  arches,  shells,  and  spherical  caps.  This  capability  represents  a 
valuable  addition  to  the  understanding  of  composite  behavior  with  material  nonlinearity. 

1.5  Current  Work 

The  subject  of  the  current  research  is  two-fold.  Previous  investigators  using  the  DSHELL 
code  have  seen  dynamic  behavior  that  is  apparently  unstable.  One  goal  of  the  current  research  is 
to  analyze  this  behavior  and  classify  its  origin  as  numerical  or  dynamical.  Some  of  the  behavior  is 
chaoslike  in  nature.  The  other  thrust  of  the  research  is  towau’d  understanding  numerical  consider¬ 
ations  in  using  the  DSHELL  program.  Key  numerical  parameters  include  the  time  step,  the 
Newton-Raphson  convergence  tolerance,  Ctoi,  and  the  coefficients  of  the  beta-m  differential  equation 
solver,  0i.  In  the  current  research,  the  only  numerical  parameter  examined  in  detail  is  the  time 
step  with  the  goal  of  devdoping  criteria  for  generating  useful  results  in  a  minimum  of  computer 
time. 

The  DSHELL  model  takes  a  considerable  amount  of  computer  time  to  run,  and  studies  of 
long-term  responses  to  loads  are  not  practical.  Because  of  this,  the  MSHELL  simplified  model  was 
developed  as  an  analog  to  the  composite  shell  model.  The  MSHELL  model  allows  for  long-term 
study  and,  by  substituting  a  sinusoidal  load  for  a  step  load,  can  mimic  the  behavior  of  a  node  in 
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the  shell’s  finite  element  model.  The  computers  used  in  the  study  were  the  SUN  SPARCstation 
330  and  the  CRAY  X-MP/216. 
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11  THEORY 


2.1  Nonlinear  Shell  Analysis 

“Finite  elements  for  shells  have  been  among  the  most  difficult  to  devise.”  (8:341)  There  have 
been  three  families  of  approaches  to  the  development  of  these  elements. 

1 .  Flat  elements:  Plane  bending  elements  with  plane  membrane  element 

2.  Curved  elements:  Classical  shell  theory 

3.  Mindlin-type  elements:  Special  (degenerated)  forms  of  solid  elements 

The  DSHELL  program,  used  in  the  current  research,  uses  an  element  which  falls  into  the  third 
category.  The  interested  reader  is  referred  to  Dennis  (9:78)  for  a  complete  development  of  the 
theory  behind  this  finite  element  code.  A  description  of  the  displacement  equations,  the  36  degree 
of  freedom  (DOF)  isoparametric  element,  and  the  numerical  solving  and  integration  schemes  are 
included  here  for  clarity. 

The  shell  geometry  is  shown  in  Figure  2.1.  The  shell  thickness  is  h  ,  the  in-plane  coordinates 
are  x  and  s  ,  and  the  depth  coordinate  z  (not  shown)  is  positive  toward  the  center  of  curvature. 
The  angle  ^  represents  the  fiber  orientation  angle,  so  that  $  =  0  represents  a  ply  whose  fibers  are 
aligned  with  the  x  (longitudinal)  aocis.  In  the  text,  subscripts  T’  and  x’  are  used  interchangeably, 
as  are  the  subscripts  ‘2’  and  ‘s'. 

2.1.1  DSHELL  Displacement  Equations.  Contracted  notation  as  described  by  Table  2.1  is 
used  in  the  following  discussion,  where  the  represent  stresses  and  the  represent  strains. 
The  DSHELL  displacement  equations  are  derived  from  the  general  curvilinear  strain-displacement 
relations  described  in  Saada  (24:136-137)  with  the  following  assumptions. 

1.  The  shell  is  in  a  modified*  state  of  plane  stress  where  try  =  0  (hence  (3  =  0). 

‘‘Modified’  in  the  sense  that  in  this  case,  through  the  thickne.ss  shear  stresses  are  not  neglected. 
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STRESS 


STRAIN 

EXPLICIT 

CONTRACTED 

fit 

ft 

(22 

(2 

(33 

^3 

2f23  =  723 


2<13 


CYLINDRICAL 

COORDINATES 


x^l 


s— 2 


z— *3 


s-z— 4 


Table  2.1.  Contracted  notation  for  stresses  and  strains 
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2.  The  transverse  shear  stresses,  (T4  and  <rs,  are  allowed  to  be  small  nonzero  values  which  vary 


parabolically  through  the  shell  thickness,  vanishing  at  the  shell  surfaces. 

3.  Only  linear  terms  are  retained  in  the  strain-displacement  equations  for  (4  and  e^. 

4.  Because  (3  =  0,  the  transverse  displacement,  w,  is  assumed  constant  through  the  shell  thick¬ 
ness. 

5.  The  ratio  of  the  shell  thickness,  h,  to  the  shell  radius  of  curvature,  R  ,  is  h/B  <  1/5. 

With  these  assumptions,  the  displacement  equations  may  be  written  as  (19:39) 


u(x,  8,  z) 

=  Uo  +  r'fl- 

-  +  "  ■ ' 

v(,x,8,z) 

=  "“['-i 

w(x,e) 

=  w 

(2.1) 


where 

•  uq  is  the  displacement  in  the  x  coordinate  direction  at  the  shell  midsurface 

•  Vo  is  the  displacement  in  the  s,  or  circumferential,  coordinate  direction  at  the  midsurface 

•  z  is  the  distance  from  the  shell  midsurface 

•  u  is  the  displacement  in  the  x  coordinate  direction 

•  V  is  the  displacement  in  the  s  coordinate  direction 

•  u)  is  the  transverse  displacement  in  the  z  direction 

•  u;,i  is  the  slope  dw/dx  at  the  node 

•  ie,2  is  the  slope  dw/ds  at  the  node 

•  is  the  rotation  of  the  normal  to  the  shell  surface  in  the  i  direction  (about  the  s  axis) 

•  ♦2  is  the  rotation  of  the  normal  to  the  shell  surface  in  the  s  direction  (about  the  x  axis) 


X 


Figure  2.2.  The  36-DOF  isoparametric  shell  element 


2.1.2  The  36~DOF  Isoparametric  Element.  The  finite  element  used  in  the  DSHELL  code 
is  derived  through  energy  methods.  The  finite  element  method  allows  the  five  coupled  nonlinear 
partial  differential  equations  representing  the  shell’s  equilibrium  to  be  reduced  to  coupled  nonlinear 
algebraic  equations  (19:69).  The  finite  element  used  in  the  DSHELL  code  is  depicted  in  Figure  2.2. 
The  isoparametric  element  has  seven  DOF  at  each  corner,  while  the  mid-side  nodes  have  only  the 
u  and  V  DOF. 

2.1.S  The  Beta-m  Method.  A  numerical  technique,  the  beta-m  method  (13)  is  employed  to 
solve  the  second-order  differential  matrix  equation 

Mx-t-Cx-t-Kx  =  F(t)  (2.2) 
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where  x  ,  x  ,  and  x  represent  the  nodal  displacement,  velocity,  and  acceleration  vectors.  The  mass 
matrix,  M  ,  is  of  the  consistent  type,  wherein  the  same  shape  functions  used  to  generate  the  element 
stiffness  matrices  are  used  to  generate  the  mass  matrix  (8:370).  The  lumped  mass  matrix  approach 
is  not  suitable  for  the  current  research,  as  rotary  inertia  would  have  to  be  arbitrarily  assigned  to 
the  nodes  to  capture  the  forces  generated  by  rotational  degrees  of  freedom.  The  damping  matrix, 
C  ,  is  generated  using  a  proportional  damping  scheme  in  which 

C  =  2^wM  (2.3) 

Where  ^  is  the  damping  ratio  and  w  is  the  linear  fundamental  natural  frequency  (6:727).  The 
displacement  vector,  x,  (for  an  r-noded  mesh),  is  of  the  form  (superscripts  denote  node  numbers, 
not  powers) 

*  =  I  u'  w\i  w^,2  ■■■  ’*'5  I 

(2.4) 

The  beta-m  method  offers  several  advantages:  (1)  it  is  a  single  step  method,  i.e.  it  needs  data 
from  only  the  most  recent  time  step  to  calculate  parameters  for  the  next  time  step,  (2)  it  is  a  very 
general  method,  i.e.  the  choice  of  parameters  0i,  /?2,  •  •  -  ,  0m-i  not  only  govern  the  stability  and 
accuracy  characteristics  of  the  method,  but  allow  the  user  to  select  from  a  broad  range  of  implicit 
or  explicit  algorithms  which  are  sub-families  to  the  beta-m  method,  and  (3)  the  finite  difference 
scheme  is  not  required  for  this  method.  The  beta-m  method  is  defined  by-  (30:383-384) 

(2.5) 

^Indices  in  parenthesis  indicate  order  of  derivative,  not  exponents. 
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where 


and 


;=* 


(2.6) 


6fc=/?th'"-‘/(m-fc)! 


(2.7) 


and  fc  =  0,  1,  . . . ,  m.  The  value  of  /?m  is  always  taken  to  be  unity,  and  h  is  the  time  increment  for 
each  time  step.  The  method  of  order  m  implies  that  xj,”**  is  the  highest  derivative  to  be  retained. 
The  terms  above  relate  the  displacement,  velocity,  and  acceleration  at  time  step  n  +  1  to  their 
values  at  time  step  n  and  the  unknown  change  in  acceleration.  Ax,  as  follows: 


Xn+l 

=  x„ -1- Xnh -1- ix„h^ -1- ^/?o/i-Ax 

(2.8) 

Xn  +  1 

=  x„ -1- Xnh -(- /?ihAx 

(2.9) 

Xn+1 

=  x„  -b  Ax 

(2.10) 

Substituting  Eq  (2.5)  into  Eq  (2.2)  at  time  tn+i  yields 


[62M  +  6iC  +  6oK(qo  +  6oAx)]  Ax  =  F„+i  -  {Mq2  +  Cqi  +  K(qo  +  6oAx)qo}  (2.11) 


where  F„+i  is  the  applied  load  vector  at  time  <n+i.  Eq  (2.11)  represents  a  set  of  nonlinear  algebraic 
equations,  and  the  Newton-Raphson  iterative  method  is  employed  to  solve  it  (12:66).  At  each  time 
step,  we  assume 

Axi+i  =  Ax,  +  ^x,  (2.12) 

where  i  is  the  iteration  number  within  a  single  time  step,  .\pplying  this  to  Eq  (2.11)  we  have 


[62M  +  61C  4-  boKr  (qo  +  6oAx,  )]  6xi  = 


F„+i  -  M  {qo  +  6oAx,}  -  C  (qi  +  61  Ax,}  -  K(qo  +  h^Ax)  {qo  +  6oAx,}  (2.13) 
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where  Kt(x)  is  a  tangentid  stiffness  matrix  evaluated  using  the  displacements  x,  which  is 


A't  (x)  =  Ao  +  A'l  (x)  +  A'2  (x) 


(2.14) 


Eq  (2.13)  is  then  solved  using  the  following  algorithm 

1.  Given  x„,  x„,  and  x„  at  time  tn,  we  seek  solutions  at  fn+i 

2.  Calculate  qo,  qi,  and  q2  from  Eq  (2.6). 

3.  Given  Ax  and  x„+i  from  the  ith  iteration,  we  obtain  the  right-hand-side  of  Eq  (2.13),  and 
the  updated  tangential  stiffness  matrix  from  Eq  (2.14). 

4.  Solve  for  the  primary  unknown,  ix,  from  Eq  (2.13). 

5.  Calculate  the  updated  solution  vector  6x,+i  from  Eq  (2.12). 

6.  Update  Xn+i,  x„+i,  and  Xn+i  for  the  (i  +  l)th  iteration  from  Eq  (2.5). 

7.  If 


+  l)i; 

.  - 

»+l 

(Xn  +  l)J^ 

. 

t 

(3Cn+l)J^ 

. 

1 

go  to  step  (1)  for  the  next  time  step,  otherwise  go  to  step  (3)  for  the  next  iteration.  The 
value  of  j  is  the  degree  of  freedom  number  and  r  is  the  total  number  of  degrees  of  freedom 
in  the  model.  The  value  of  c  is  chosen  by  the  user  and  represents  the  maximum  acceptable 
convergence  tolerance. 


2.8  A  Simplified  Shell  Model 

The  DSHELL  finite  element  model  incorporates  many  complex  nonlinear  phenomena  involv¬ 
ing  hundreds  or  even  thousands  of  active  and  interacting  degrees  of  freedom.  Over  ten  thousand 
lines  of  code  are  executed  to  perform  just  one  iteration  in  one  time  step,  using  large  amounts  of 
computer  memory  and  time.  It  would  be  desirable  to  have  a  simpler  analogy  to  use  cis  a  tool  to 
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understand  some  of  the  results  produced  by  DSHELL.  To  this  end,  another  computer  program  was 
developed  to  try  to  model  some  of  the  important  features  of  shell  behavior  on  a  smaller  and  more 
easily  analyzed  scale. 

The  model  described  herein,  referred  to  as  ‘MSHELL’,  incorporates  some  of  the  important 
features  of  the  laminated  cylindrical  shell: 

•  The  ‘snap-through’  phenomena  is  facilitated. 

•  Stiffness  exists  in  both  in-plane  directions,  as  well  as  in  the  transverse  direction. 

•  The  equations  of  motion  for  this  system  are  nonlinear  though,  as  in  the  finite  element  program, 
the  material  properties  (in  this  case,  the  spring  stiffness)  are  linear. 

•  The  coupled  equations  meet  the  criteria  for  admitting  chaotic  behavior. 

The  model  does  not  mimic  the  cylindrical  shell  in  these  important  ways: 

•  The  differences  in  material  properties  through  the  laminate  thickness  are  not  modeled,  as  no 
thickness  related  parameters  exist 

•  The  shell’s  resistance  to  transverse  displacement  due  to  its  width  (in  the  longitudinal  direc¬ 
tion)  is  not  modeled.  (The  model  is  more  ‘arch-like’  than  ‘shell-like’  in  this  regard.) 

•  The  distributed  mass  of  the  cylindrical  shell  is  not  modeled. 

•  The  interactions  of  the  motions  at  various  points  is  not  modeled,  though  interactions  of  nodes 
in  the  suddenly  loaded  shell  may  be  effectively  simulated  through  sinusoidal  loading  of  the 
model. 

2.S.1  Equations  of  Motion.  The  equations  of  motion  are  purposefully  derived  from  kine¬ 
matic  considerations  rather  than  energy  because  it  is  desired  to  use  energy  balance  as  a  checking 
mechanism. 


2-8 


Figure  2.3.  Simplified  model  of  shell 

The  model  is  depicted  in  Figure  2.3,  and  its  appearance  as  it  relates  to  the  shell  geometry  is 
shown  in  Figure  2.4.  All  springs  are  of  undeformed  length  6,  and  springs  1  and  2  remain  horizontal 
(in  the  x  —  y  plane)  regardless  of  the  deformed  geometry,  such  that  the  ends  of  springs  1  and  2  can 
be  thought  of  as  riding  in  tracks  of  infinite  length  offering  no  resistance  to  vertical  motion  of  the 
spring  ends.  The  goal  here  was  to  try  to  keep  the  model  as  simple  as  possible  while  still  capturing 
the  nonlinear  effects.  Attempting  to  introduce  additional  transverse  stiffness  due  to  springs  1  and 
2  severely  complicated  the. model,  though  it  should  be  noted  that  this  feature  explicitly  decouples 
longitudinal  and  transverse  stiffnesses.  The  object  is  of  mass  m  ,  and  the  dimensions  h  and  d  are 
the  projections  of  undeformed  spring  3  (or  4)  on  the  horizontal  and  vertical,  respectively. 
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Figure  2.4.  Simplified  model  of  shell  as  it  relates  to  geometry  of  shell,  radius  r 


The  equations  of  motion  for  the  model  are 


u(u,v,w,u,fu{t))  = 


(1/^)  j/u  +  ^*1 - j - l>^2  + 

,  /  1  I  ^^1  bki  bks  bk^  \ 

+  (  - ^2  +  "7 - fcs  +  ‘T - ^4  +  “7“  I  U  —  CU 

\  0l  *2  03  04  / 

t;(u,r,u;,  V, /„(<))  = 

(l/ra)  [/.  -  Uj  +  *^  +  W,  - 

V  61  *2  63  64  / 

w(u,v,w,w,U{t))  = 


(1/m) 


f  ,<1.  ,  ,,  Wi4  ^^3  ,  bk4\ 

fw  +  dkz - 7 - h  dk^ - r - h  (  — A:3  +  - k4  +  - —  j  li?  —  cit; 

h  04  \  63  64  / 


(2.16) 


(2.17) 


(2.18) 


where  the  variables  are  defined  as  follows; 
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•  u,  V,  and  w  are  the  displacements  in  the  x,  y,  and  z  coordinate  directions  respectively. 

•  The  dots  represent  derivatives  with  respect  to  time,  e.g.  u  =  dufdt. 

•  The  ki  represent  the  linear  spring  stiffnesses. 

•  6  is  the  undeformed  spring  length,  and  is  the  same  for  all  springs. 

•  The  bi  are  the  deformed  spring  lengths,  the  equations  for  which  are  found  in  Eq  (A.l)  on 
page  A-1. 

•  ft  is  the  projection  of  undeformed  spring  3  (or  4)  on  to  the  y  coordinate  direction. 

•  d  is  the  projection  onto  the  z  coordinate  direction. 

•  c  is  the  viscous  damping  coefRcient 

•  The  /,  are  the  time  varying  components  of  the  applied  force. 

The  derivation  of  the  equations  of  motion  is  provided  as  Appendix  A . 


2.2.2  Static  Collapse  Load.  Calculation  of  the  static  collapse  load  is  not  difficult  with  this 
model.  The  total  static  spring  force  in  the  w  direction  due  to  springs  3  and  4  (in  the  presence  of 
only  a  w  displacement)  is  given  by 


-2k 


-b+yjh^-d-id-wf 


(d  —  w) 


/spring 


+  (d-  w)^ 


(2.19) 


Its  derivative,  found  using  the  D  command  in  Mathemattca^  (3:52),  is 


_d 

dw 


—  f 


2k  (-b+  yjh^-  +  {d-w)-'^ 


2k 


ijh^  +  (d-  wf 

^-6+v^ 


+  (d  -  w)  (d  -  w) 


(^h^-+(d-w)-y 


+ 


^  Mathematica  is  a  registered  trademark  of  Wolfram  Research.  Inc. 
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2k{d-wf 
’h^  +  (d-wf 


The  real  roots  of  this  derivative,  using  the  Solve  command  in  Mathematica  (3:46)  are 


lUl  = 


2<i  +  2\/6^ 


W2  = 


2d-2\/bi  -  hih^ 


(2.20) 


The  displacement  at  the  static  snapping  load  is  given  by  W2,  for  exceeding  this  displacement  moves 
us  to  the  unstable  portion  of  the  equilibrium  curve  (dP/dw  <  0): 


.  2d-2\/bi  -h^hi 

w{fmax)  =  - :: - 


(2.21) 


Substituting  this  value  into  Eq  (2.19)  yields 


2\/bi-hihi  (6  -  \/bhTj  k 


(2.22) 


where  fmap  represents  the  static  collapse  load  in  the  presence  of  only  transverse  (no  in-plane) 
loading.  The  static  collapse  load  is  seen  to  be  a  function  of  only  the  undeformed  spring  length,  the 
‘shell’  depth,  and  the  spring  stiffness.  Eq  (2.19)  is  plotted  in  Figure  3.1  on  page  3-2. 


S.S.S  Dynamic  Collapse  Load.  Calculation  of  the  dynamic  collapse  load  must  include  two 
quantities  ignored  in  the  calculation  of  the  static  collapse  load:  mass  (inertia)  and  the  load  function 
/(t)  .  Eq  (2.18)  is  the  differential  equation  representing  the  motion  in  the  w  direction.  In  an  attempt 
to  find  a  closed  form  solution,  we  restrict  the  motion  to  the  w  axis,  i.e.  u  =  v  =  0,  and  we  exclude 
damping.  Under  these  conditions,  Eq  (2.18)  becomes 

mii)  -f-  2kw  +  k  i^b{d  -  w)  “  2d|  =  /„,  (2.23) 
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But  when  u  =  i>  =  0, 63  =  64  and  Eq  (2.23)  becomes 


mw  +  2kw  +  k  |6(d  -  w)  ~  2d|  =  /u, 


which  simplifies  to 


mw 


+  2k(cl 


Eq  (A.l)  gives  us  63  which,  with  u  =  t;  =  0  is 


63  =  \//»^  +  (d—  wY 


Substitution  yields 

+  (2.24) 

Despite  the  simplification  of  confining  the  motion  to  the  z  axis,  Eq  (2.24)  is  a  non-linear  second 
order  differential  equation  whose  closed  form  solution  is  not  obvious  and  likely  exists  only  for  certain 
parameters.  Mathematical  theorems,  such  as  the  uniqueness  and  existence  theorem  (5:24-25)  reveal 
the  conditions  for  which  solutions  are  guaranteed  to  exist,  but  do  not  produce  those  solutions.  For 
describing  the  general  motion  of  the  mass,  numerical  solutions  may  be  had  more  quickly  than 
analytical  ones.  Furthermore,  under  cyclic  loading,  the  dynamic  collapse  load  loses  its  meaning  in 
that  there  is  no  ‘threshold  value’  at  which  it  may  be  said  instability  (collapse)  is  sure  to  occur  (see 
Section  3.4,  page  3-11). 

2.S.4  A  Checking  Mechanism.  It  would  be  desirable  to  have  an  ‘independent’  check  that 
the  computer  program  is  converging  to  a  reasonable  solution  for  each  time  increment  At.  Toward 
this  goal,  a  checking  mechanism  is  incorporated  into  the  model  to  monitor  the  energy-balance  for 
the  no  damping  case.  The  idea  for  using  energy  as  a  check  comes  from  Kounadis’  work  in  studying 
chaoslike  phenomena  of  suddenly  load  structures  (14:302). 
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The  check  is  performed  as  follows: 


1.  The  cumulative  error,  E  ,  the  integrated  work  performed,  W  ,  and  thj  cumulative  percent 
error,  Epct  are  set  equal  to  zero: 


E  =  W  =  Ep^t  =  0 


(2.25) 


2.  At  the  end  of  the  ith  time  step,  the  change  in  kinetic  energy,  6Ti  ,  potential  energy,  SV]  , 
and  work  performed  by  the  external  force,  8Wi  (using  the  average  force  of  the  current  and 
previous  time  step)  are  calculated: 


STi 


8Vi 


8Wi 


Lf=i 

3 


(=1 


J=1 


(2.26) 

(2.27) 

(2.28) 


where  the  fcj  are  the  spring  stiffnesses  for  each  of  the  four  linear  springs,  the  6j  are  the 
deformed  spring  lengths,  b  is  the  undeformed  spring  length  (equal  for  all  springs),  m  is  the 
mass,  and  the  j  index  refers  to  the  component  directions,  e.g.  £2  =  v  and  fa  =  fw 

3.  Next,  the  error  for  the  ith  time  step,  8Ei  ,  and  is  calculated: 


8Ei  =  \6Wi  -  6Ti  -  6Vi\ 


(2.29) 


V 


4.  The  cumulative  error  and  cumulative  work  performed  may  then  be  calculated,  yielding  the 
cumulative  percent  error: 


Ei  —  Ei- 1  +  8Ei 
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Wi  =  \V:.l+6Wi 

(£,«),  =  >00  §■ 

(2.30) 

calculated  as  follows: 

/.  \  _  .jYln=l^^Pot)n 

(ermsU  -  y  j 

(2.31) 

At  the  end  of  each  time  step,  Crm»  is  compared  to  a  user  defined  meiximum  (typical  values  are 
employed  in  Section  3.4).  If  the  maximum  error  is  exceeded,  the  program  notifies  the  user  and 
terminates,  otherwise,  another  lime  step  is  taken  and  the  error  checked  beginning  with  step  2.  The 
FORTRAN  implementation  of  the  energy  balance  scheme  is  provided  in  Appendix  C. 

S.S.5  The  Integration  Scheme.  (31:190-193) 

An  iterative  procedure  is  used  to  solve  the  nonlinear  equations  of  motion.  This  predictor- 
corrector  method  involves  predicting  the  outcome  for  a  step  in  time,  then  applying  an  implicit 
correcting  formula  to  improve  the  result.  The  subsequent  result  is  then  used  as  the  prediction,  the 
corrector  is  applied,  and  the  process  continues.  This  particular  algorithm  was  chosen  for  several 
reasons:  (1)  the  method  is  specifically  designed  for  second-order  ordinary  differential  equations,  (2) 
it  is  an  implicit  amd  unconditionally  stable  algorithm,  and  (3)  with  only  three  degrees  of  freedom, 
employing  a  more  complex  matrix-oriented  scheme  is  not  necessary  or  practical  —  the  three  closed- 
form  solutions  for  the  accelerations  are  easily  calculated  by  the  computer,  eliminating  the  need  for 
iterating  to  solve  for  accelerations  (as  in  the  beta-m  method).  The  algorithm  will  be  described  for 
the  w  coordinate,  though  calculations  for  the  u  and  v  coordinates  are  identical. 

For  the  first  iteration,  the  velocity  is  estimated  at  the  end  of  the  first  time  step: 


tei  =  leo  +  woAto 


(2.32) 
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Where  wq  and  vbo  are  the  initial  velocity  and  acceleration  in  the  w  direction.  The  displacement  w 
is  then  predicted  based  upon  the  predicted  velocity; 


wi  =  Rw  +  -ti;iAto 


(2.33) 


Where  Ry,  =  ivq  +  ^wq^I,  with  wq  being  the  initial  displacement.  The  closed  form  solution  for  the 
acceleration,  shown  in  Eq  (2.18),  is  then  used  to  calculate  w. 

For  subsequent  iterations  i  on  time  step  j,  j  >  1,  the  procedure  is  similar.  First,  velocity  is 
predicted: 

=  {Qw)j-i  +  ^(u)j)i-iA<j_i  (2.34) 

Where  (<3ui);-i  =  Then,  displacement  is  predicted  based  upon  the  predicted 

velocity: 

(u/j)j  =  (Rv/)j-i  +  (2.35) 


Convergence  is  checked  for  by  comparing  the  difference  between  the  displacement  calculated  for 
the  current  iteration  and  the  displacement  calculated  in  the  previous  iteration.  If  the  difference  is 
less  than  some  arbitrary  user-defined  c,  the  iteration  is  successful.  If  not.  a  new  iii  is  calculated 
and  a  new  iteration  is  begun.  As  in  the  DSHELL  finite  element  model,  MSHELL  uses  an  averaging 
scheme  to  check  for  convergence,  such  that,  at  the  ith  iteration. 


>)’) 

1-1 

<_e 

i 

(2.36) 


must  be  satisfied  for  convergence  to  be  deemed  successful  (30:2-38). 

Like  the  beta-m  method,  this  is  a  single-step  method,  allowing  change  of  the  time  step  during 
a  computer  run.  The  FORTRAN  implementation  of  this  method  is  provided  as  Appendix  B. 
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Chaos 


e.s 

There  Me  three  classic  types  of  dynamical  motion: 

1.  Equilibrium 

2.  Periodic  motion  or  a  limit  cycle 

3.  Quasiperiodic  motion 

In  phase  space,  the  equilibrium  state  is  associated  with  a  point,  the  periodic  motion  with  a  closed 
curve,  and  the  quasiperiodic  motion  with  a  surface  in  'hree-dimensional  phase  space.  These  states 
are  called  attractors,  since,  with  damping,  the  dynamical  .■;> stem's  motion  tends  to  decay  to  one 
of  these  three  states.  Chaotic  behavior,  however,  involves  yet  another  attractor  called  a  strange 
attractor.  Unlike  the  classic  attractors,  this  attractor  is  not  associated  with  a  classical  geometric 
object  in  phase  space,  but  a  new  geometric  object  called  the  fractal  set.  In  three-dimensional  phase 
space,  the  fractal  set  of  a  strange  attractor  looks  like  “. . .  a  collection  of  an  infinite  set  of  sheets  or 
parallel  surfaces,  some  of  which  are  separateu  by  distances  that  approach  the  infinitesimal  ( 17:23)  " 

There  are  several  useful  graphical  tools  for  the  analysis  of  nonlinear  systems.  Consider  a 
linear^,  one  degree  of  freedom,  forced  spring-m^lss  system  (31:48-49).  The  motion  includes  the 
steady  state  response  at  the  frequency  of  the  forcing  function,  Q  ,  superimposed  on  the  transient 
response  at  the  natural  frequency  of  the  system,  w  =sfkfm.  The  second  order  differential  equation 
of  motion  for  this  system  is 

mx  +  kx  =  f  (2  37) 


and  its  solution  for  z(0)  =  i(0)  =  0  is 


(2.3S) 


*  A  linear  system  is  used  here  for  illustration  as  the  closed  form  solution  of  the  differential  equation  is  immediately 
available. 
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Figure  2.5.  Displacement  vs.  time  for  forced,  undamped  linear  vibration 

The  motion  versus  time  for  a  particular  choice  of  k,  m,  f,  and  f2  is  depicted  in  Figure  2.5.  The 
phase  diagram  plots  velocity  versus  displacement  for  a  dynamical  degree  of  freedom.  The  path 
traced  by  the  function  is  referred  to  as  the  trajectory.  The  phase  diagram  for  the  example  system 
is  shown  in  Figure  2.6. 

The  Poincare  section  plots  points  which  represent  snapshots  in  time  of  the  phase  diagram. 
The  phase  diagrcun  is  sampled  at  a  particular  frequency  (usually  the  frequency  of  the  forcing 
function)  amd  the  sampled  points  form  the  Poincare  section.  The  Poincare  map  is  also  referred  to 
as  a  surface  of  section,  because  it  may  be  generated  by  sampling  the  displacement  and  velocity  as 
the  object  crosses  a  plane  or  surface  in  space.  Among  other  things,  the  Poincare  section  indicates 
resonances  between  the  forcing  function  and  the  dynamical  system’s  response.  The  Poincare  map 
for  the  example  system  is  shown  in  Figure  2.7.  The  sampling  rate  is  the  frequency  of  the  forcing 
function,  with  a  phase  angle  of  zero. 
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Figure  2.6.  Phase  diagram  for  forced,  undamped  linear  vibration 
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Figure  2.7.  Poincare  section  for  forced,  undamped  linear  vibration 
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The  Fourier  transform  may  be  used  to  determine  the  power  spectrum  of  a  particular  motion. 
The  power  spectrum  relates  frequencies  of  motion  to  the  energy  the  frequencies  contain.  Power 
spectra  with  discrete  peaks  of  energy  indicate  non-chaotic  behavior,  while  broad-banded  ‘noise-like’ 
spectra  may  indicate  chaotic  behavior.  Near  the  onset  of  chaos,  the  discrete  peaks  are  absorbed 
into  a  more  continuous  distribution  of  frequencies.  With  fully  chaotic  behavior,  the  continuous 
spectrum  may  dominate  the  discrete  spikes  (17:148).  The  Fast  Fourier  Transform  (FFT)  is  used  in 
this  research.  It  has  the  advantage  of  needing  on  the  order  of  N  logj  N  operations  as  compared  to 
A’^  operations  for  the  Discrete  Fourier  Transform  (DFT),  where  N  is  the  number  of  sampled  data 
points  (11:766).  In  the  DFT  (and  hence  the  FFT,  which  is  but  a  special  algorithm  for  calculating 
the  DFT)  the  sampled  data  set  is  represented  by  a  set  of  Fourier  coefficients  Ar{u)) 

N-l 

"^'■(w)  =  ^  Xt  exp{-2irirk/N),  r  =  0, 1 . A"  -  1  (2.39) 

k=0 

where  X*  is  the  sampled  data  set  (in  this  case,  the  displacement  samples),  and  i  =  y/^.  The 
coefficients  >lr(w)  are  generally  complex,  so  it  is  useful  to  define  a  real- valued  function  called  the 
power  spectrum,  5(w)  =((ylr(w)|p  (2:31).  The  FFT  is  used  to  generate  Figure  2.8,  which  reveals  the 
forcing  frequency,  /n  =  250  Hz,  and  the  system’s  natural  frequency,  fu  =  1597  Hz.  The  amplitude 
information  reveds  how  relatively  little  excitation  is  caused  when  the  forcing  frequency  is  far  from 
the  system’s  natural  frequency.  A  key  rule  of  thumb  in  taking  data  from  the  time  domain  to  the 
frequency  domain  is  that  the  sampling  rate  should  be  at  least  twice  that  of  the  highest  frequency 
of  interest.  In  this  research,  the  time  steps  used  are  less  than  0.0001  seconds,  corresponding  to 
a  sampling  rate  of  10  kHz,  hence  capturing  frequencies  below  5  kHz.  And  for  many  plots,  the 
sampling  rate  is  much  higher. 

The  bifurcation  diagram  plots  velocity  versus  some  other  system  parameter  (such  as  the  am¬ 
plitude  of  the  applied  load).  The  points  are  sampled  as  in  the  Poincare  section.  Hence  each  vertical 
slice  of  the  bifurcation  diagram  may  be  thought  of  as  a  collapsed  Poincare  section  for  a  particular 
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Frequency  (Hz) 


Figure  2.8.  Relative  power  vs.  frequency  for  forced,  undamped  linear  vibration 


load  condition.  The  bifurcation  diagram  depicts  the  solution(s)  to  the  nonlinear  differential  equa¬ 
tion  given  a  particular  forcing  function.  In  generating  the  bifurcation  diagrams,  care  must  be  taken 
to  discard  early  transient  motion  (17:130).  Figure  3.12  on  page  3-10  is  a  bifurcation  diagram. 

Two  specific  types  of  bifurcation  are  seen  in  the  current  research.  Period  doubling,  also  called 
the  flip  bifurcation  occurs  when  different  oscillatory  motions  of  the  same  frequency,  w,  combine  in 
such  a  fashion  as  to  repeat  at  a  frequency  of  w/2  (2:72).  When  this  happens,  the  closed  orbit  (in 
phase  space)  representing  the  fundamental  frequency  loses  its  stability  and  is  replaced  by  another 
(stable)  closed  orbit  whose  period  is  twice  that  of  the  original  (29:127).  The  flip  bifurcation  is 
peculiar  to  cycles  and,  unlike  the  fold  bifurcation,  has  no  correspondence  for  equilibria  (29:127). 

The  fold  bifurcation,  also  called  a  saddle-node  bifurcation,  is  exemplified  in  the  snap-buckling 
of  a  support  arch  (or  shell)  under  vertical  load  (29:256).  The  fold  bifurcation  is  a  structurally  stable 
bifurcation,  arising  when  a  pair  of  equilibria  exist  —  one  unstable  and  the  other  attracting  (29:256). 
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2-4  Selection  of  Time  Step,  At 


Choosing  m  appropriate  time  step,  At  ,  for  a  given  dynamic  nonlinear  problem  is  both  very 
important  and  very  difficult.  Both  the  behavior  of  the  forcing  function,  f(t),  and  the  behavior 
(response)  of  the  dynamical  system  acted  upon  by  f(t)  must  be  captured  (15:550).  Depending  on 
the  integration  technique,  the  choice  of  At  may  determine  the  stability  of  the  integrating  scheme. 
Even  in  those  cases  where  unconditionally  stable  integrating  techniques  are  used  (such  as  in  the 
present  research),  numerical  accuracy  is  not  assured — the  algorithm  may  always  converge,  but  it 
may  converge  to  a  solution  not  of  the  desired  accuracy.  Unfortunately,  decreasing  the  time  step 
does  not  guarantee  an  improvement  in  accuracy.  There  is  a  minimum  At  below  which  accuracy 
is  compromised  due  to  increasing  the  number  of  calculations  performed  over  the  duration  of  the 
computer  run  (round-off  error).  In  addition,  practical  considerations,  such  as  limits  on  computer 
time,  drive  the  analyst  toward  using  the  largest  possible  time  step  that  gives  acceptable  results. 
The  DSHELL  program  executes  over  10,000  lines  of  code  for  a  single  iteration  of  a  single  time 
step  (27),  and  it  is  not  uncommon  for  5  or  more  iterations  to  be  required  for  convergence  in  a  time 
step.  For  most  of  the  runs  done  in  the  present  research,  an  extremely  small  convergence  tolerance 
of  (10~^°)  percent  was  used.  This  was  done  in  an  effort  to  eliminate  cumulative  error  due  to  poor 
accuracy  in  the  convergence  as  a  source  of  instability.  Happily,  the  choice  of  such  a  small  tolerance 
does  not  significantly  increase  the  number  of  iterations  required  for  convergence  —  the  choice  of 
At  has  a  much  more  pronounced  effect. 

Many  investigators  have  suggested  various  ‘rules  of  thumb’  for  choosing  an  appropriate  At 
Table  2.2  lists  some  of  these  rules  and  their  sources.  In  Table  2.2,  ri  is  the  period  associated 
with  the  first  fundamental  frequency  of  vibration,  Tmax  is  the  period  associated  with  the  highest 
frequency  of  interest,  emd  r  represents  the  frequency  of  oscillation  of  a  single  degree  of  freedom 
system.  In  addition  to  these,  Almroth  and  Brogan  (1:6-32-6-37)  offer  several  methods  with  the 
caveat  that  trial  and  error  may  be  required  following  the  first  choice  for  At.  Almroth,  Brogan. 
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1  Technique 

Source 

Low  (15:556) 

^^max  —  ^max/20 

or  less  for  nonlinear  problems 

Cook,  et  al.  (8:407) 

Katona  and  Zienkiewicz  (13:1356) 

Table  2.2.  Techniques  for  choosing  At 


and  Cook  note  well  that  modification  of  the  initial  choice  for  At  may  be  required,  particularly  in 
the  nonlinear  case.  Employing  any  of  these  'rules  of  thumb’  for  choosing  an  appropriate  time  step 
involves  having  a  priori  knowledge  of  how  the  system  behaves. 

What  makes  the  choice  of  time  step  particularly  difficult  is  the  i)hysical  behavior  of  the 
nonlinear  problem.  For  example,  the  inertial  forces  during  and  after  collapse  (snapping)  may  be 
much  larger  than  before  collapse  or  during  steady-state  oscillation.  An  eigenvalue  analysis  yields 
clues  as  to  the  frequencies  of  interest,  but  only  in  the  neighborhood  of  the  undeformed  geometry. 
The  eigenvalues  (and  hence  the  natural  frequencies)  change  as  the  object  deforms.  In  the  present 
research,  a  At  suitable  for  prebuckling  or  steady-state  response  was  unsuitable  for  postbuckling 
behavior.  In  studies  of  the  laminated  composite  arch,  the  postbuckled  deformed  geometry  was 
seen  to  change  radically  over  just  one  step  in  time,  casting  doubt  on  whether  the  geometry  found 
was  correct  or  just  another  (incorrect)  geometry  satisfying  the  energy  balance  equations  of  the 
structure. 

2.5  Post-processing  Capability 

As  part  of  this  research,  two  FORTRAN  routines  were  written  to  facilitate  post-processing 
of  the  DSHELL  output. 

Since  Dennis  (9)  authored  the  SHELL  code  in  1988,  users  of  the  code  have  used  various 
means  of  analyzing  the  results.  However,  no  full-featured  post-processing  capability  was  available. 
Minor  changes  have  been  made  to  the  code  at  various  times  to  allow  plotting  programs,  such 
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as  GRAPHER™,  to  present  limited  amounts  of  data  in  a  more  easily  evaluated  form,  but  no 
means  have  been  available  to  view  the  entire  deformed  structure  and  the  associated  displacements, 
velocities,  accelerations,  strains,  and  stresses  in  a  comprehensive  way.  For  this  reason,  one  goal  of 
this  research  was  to  incorporate  such  a  capability. 

2.5.1  Subroutine  UNV.  The  FORTRAN  routine  UNV  allows  the  user  of  DSHELL  to  capture 
most  of  the  output  generated  by  the  program.  UNV  generates  a  file®  in  a  format  SDRC’s  I-DEAS® 
computer  program  can  read.  The  extensive  post-processing  capabilities  of  I-DEAS  allow  the  user 
to  view  the  data  generated  by  DSHELL  in  a  graphical  way.  The  user  may  view  the  structure  in 
perspective  views  and  see  the  actual  deformed  geometry.  The  user  may  view  contour  plots  (in  color) 
of  various  quantities  (velocity,  stress,  rotational  acceleration,  etc.)  superimposed  on  the  deformed 
geometry. 

I-DEAS  also  gives  the  user  the  ability  to  plot  the  data  from  the  DSHELL  run  and,  if  the 
plotting  style  of  another  plotter  is  preferred,  I-DEAS  can  output  the  data  in  tabular  (spreadsheet) 
form  for  use  by  another  plotting  program  of  the  user’s  choice. 

The  user  of  UNV  may  specify  any  or  all  of  the  following  data  to  be  captured: 

•  nodal  displacements  (including  rotations) 

•  nodal  velocities 

•  nodal  accelerations 

•  elemental  stresses 

•  elemental  strains 

If  the  particular  problem  can  be  represented  using  one-fourth  of  a  bilaterally  symmetric 
structure,  the  user  may  choose  the  symmetry  option  when  running  UNV.  Choosing  this  option 

®SDRC  refers  to  this  as  a  ‘universal’  file. 

*I-DEAS  and  SDRC  are  trademarks  of  Structural  Dynamics  Research  Cornoration. 
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causes  UNV  to  reflect  the  data  from  the  article  into  three  adjacent  quadrants,  allowing  the  user 
to  ‘see’  a  representation  of  the  entire  article.  This  option  is  only  successful  when  Node  Number  J 
resides  at  the  intersection  of  the  two  axes  of  symmetry. 

UNV  automatically  chooses  either  rectangular  Cartesian  or  cylindrical  coordinates  depending 
on  the  model  chosen.  If  an  arch  or  shell  is  modeled,  UNV  will  generate  'lata  in  cylindrical  coordi¬ 
nates.  Choosing  a  plate  model  causes  UNV  to  generate  results  in  rectangular  Cartesian  coordinates. 
In  any  case,  the  graphics  presented  by  I-DEAS  will  depict  a  rectangular  system  in  the  corner  of 
the  terminal  screen.  For  cylindrical  coordinates,  the  axis  directions  on  the  screen  correspond  to  the 
article’s  coordinate  system  as  follows: 

•  the  +x  direction  corresponds  to  +r  {—w) 

•  the  +y  direction  corresponds  to  +0  (-fs) 

•  the  +2  direction  corresponds  to  -t-2  (+u) 

In  DSHELL,  w  displacements  for  the  mid-side  nodes  are  implied,  not  calculated.  In  order 
to  plot  the  displacements,  UNV  generates  a  cubic  curve  fit  that  matches  the  slope  {w,i  or  ii;  2) 
at  the  corner  nodes  and  has  the  appropriate  calculated  u  or  v  displacement  of  the  mid-side  node 
(depending  on  the  direction  of  interpolation). 

2.5.2  Subroutine  UMODAL.  The  FORTRAN  routine  UMODAL,  allows  the  user  to  view, 
using  I-DEAS,  the  physical  mode  shapes  associated  with  the  eigenvalues  calculated  by  DSHELL. 
DSHELL  calculates  the  displacements  and  rotations  at  each  element  corner  node,  but  UMODAL 
captures  only  the  displacements,  not  the  rotations. 

2.5.3  The  Input  Deck.  Use  of  the  post-processing  features  requires  inserting  one  card  into 
the  input  deck.  It  becomes  the  first  card,  and  is  of  the  form 

MAKE.UIV .  ISYMM .  U.LIMIT ,  MB.MAX ,  MO.PRIIT ,  U.DISP .  U.VEL ,  U_  ACC ,  U.STRESS ,  U.STRAII 
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where  the  values  of  the  variables  and  their  functions  are  described  in  Table  2.3. 


Variable 

Data  Type 

Function 

NAKEJIIV 

infeyer 

=  0,  do  not  create  universal  file 
=  1,  create  universal  file 

ISYMM 

integer 

=  0,  do  not  use  symmetry 

=  1,  use  symmetry  by  reflecting  the  model  about  Node  1  assuming 
bilateral  symmetry 

UXIMIT 

integer 

=  0,  do  not  limit  the  size  of  the  universal  file 
=  1,  limit  the  size  of  the  universal  file  to  MB.MAX  megabytes 

NBJMX 

integer 

(see  U XIMIT) 

NO-PRIIT 

integer 

=  0,  print  stress  data  to  output  file'^“"^ 

=  1,  do  not  print  stress  data  to  output  file 

UJ)ISP 

integer 

=  0,  do  not  capture  nodal  displacements  and  rotations 
=  1,  capture  nodal  displacements  and  rotations 

U-VEL 

integer 

=  0,  do  not  capture  nodal  velocities 
=  1,  capture  nodal  velocities 

UJICC 

integer 

=  0,  do  not  capture  nodal  accelerations 
=  1,  capture  nodal  accelerations 

U^TRESS 

integer 

=  0,  do  not  capture  elemental  stresses 
=  1,  capture  elemental  stresses”^®*^® 

U^TRAII 

integer 

=  0,  do  not  capture  elemental  strains 
=  1,  capture  elemented  strains 

Note:  If  the  value  of  (UTKESS  is  1 ,  then  the  user  must  select  ISTKES  to  equal  the  total  number  of  elements 
and  ISTRES(I)  must  contain  each  and  every  element  number.  This  will  generate  very  large  output  files  unless 
■0  J>RIIT  is  set  to  one. 

Table  2.3.  Variable  list  for  data  input  card  for  use  with  postprocessor 
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Ill  MSHELL;  THE  SIMPLIFIED  SHELL  MODEL 


The  MSHELL  model  was  generated  with  the  goal  of  getting  insight  into  some  of  the  phe¬ 
nomena  observed  through  the  DSHELL  finite  element  code.  The  MSHELL  code  is  only  a  small 
fraction  of  the  size  of  DSHELL,  and  run  times  are  commensurately  smaller.  Through  the  use  of 
MSHELL,  thousands  of  oscillations  may  be  observed  in  a  few  minutes  of  computer  time,  whereas 
with  DSHELL,  ten  oscillations  of  the  vibrating  shell  may  take  several  days  of  computer  time  to 
generate.  The  MSHELL  model  by  no  means  has  the  fidelity  of  the  finite  element  model,  but  many 
of  the  fundamental  characteristics  of  shell  behavior  are  modeled  (see  page  2-8). 

S.l  The  Model  Parameters 

Parameters  for  the  model  are  chosen  to  match  the  geometry,  linear  fundamental  natural 
fre<|uency,  and  static  buckling  load  of  the  cross-ply  laminated  shell  used  in  the  present  research. 
The  stiffness  of  springs  3  and  4  determine  the  static  collapse  load  (with  w  displacement  alone)  of 
the  MSHELL  model.  The  constants  of  springs  1  and  2  were  empirically  chosen  to  match  the  ratio 
of  in-plane  to  transverse  amplitudes  seen  in  edge  nodes  of  the  composite  shell  subjected  to  a  sudden 
load.  The  parameters  are  listed  in  Table  3.1.  These  parameters  yield  the  linear  natural  frequencies 


ki,  Ibs/in 

2 

3 

4 

b,  in 

h,  in 

d,  in 

m,  slugs 

39,235 

39,235 

78,470 

78,470 

5.9377 

5.7532 

1.4686 

1.5588  (lO--*) 

Table  3.1.  Parameters  for  simple  shell-like  model 


shown  in  Table  3.2.  The  first  mode  is  seen  to  be  the  transverse  vibration,  and  the  second  the 
in-plane  vibration.  The  load  versus  displacement  described  by  Eq  (2.19)  for  this  model  is  depicted 
in  Figure  3.1. 
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Table  3.2.  Linearized  behavior  of  simple  shell-like  model 


Figure  3.1.  Load  vs.  displacement  for  simple  shell-like  model.  The  static  snap  load  is  depicted  at 
point  ‘A’. 
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Figure  3.2.  Bifurcation  plot  for  MSHELL  model,  sudden  loading,  r  =  10®  Ibs/sec,  depicting 
collapse  load  at  A 


S.2  Response  to  Sudden  Loading 

The  response  of  the  MSHELL  model  to  sudden  loading  was  observed  at  two  different  load 
rates,  r  ,  both  of  which  result  in  collapse  (buckling)  of  the  model.  The  load  is  applied  in  a  highly 
eccentric  fashion  in  order  to  excite  in-plane  as  well  transverse  motion.  Figure  3.2  is  the  bifurcation 
diagram  for  the  suddenly  loaded  model.  This  plot  is  generated  by  forming  the  Poinc^l^e  map  at 
each  amplitude  and  collapsing  it  in  the  displacement  direction.  In  other  words,  each  vertical  slice 
of  the  bifurcation  diagram  is  a  Poincare  map  collapsed  in  the  z  direction.  The  arrow  at  point  A 
indicates  the  onset  of  collapse  behavior. 


Figure  3.3.  MSHELL  model  w  displacement  vs.  time  for  step  load  of  46,000  lbs  applied  at  10® 
Ibs/sec 


Id  the  first  case,  the  model  is  loaded  at  a  rate^  of  r  =  10®  Ibs/sec: 

{10®f,  <<  46,000/10® 

(3.1) 

46,000,  <>  46,000/10® 

1 

Here  f  is  the  load  function.  Under  these  conditions  the  model  collapses  eind  assumes  a  steady  stati- 
oscillation  with  w  >  d.  The  displacement  versus  time  plot  may  be  seen  in  Figure  3.3.  In  phase 
space  (Figure  3.4),  this  undamped  system  is  seen  to  rapidly  converge  to  a  fixed  orbit  —  a  purely 
periodic  motion.  The  frequency  response  generated  by  these  conditions  is  depicted  in  Figure  3.5, 
with  the  energy  peak  at  228  Hz.  The  subsequent  peaks  represent  the  second  and  tliird  harmonics 
of  this  fundamental  nonlinear  frequency  of  vibration. 

‘This  load  rate  is  chosen  because  it  is  frefinently  used  in  the  finite  element  model  analyses. 
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Figure  3.5.  MSHELL  model  frequency  response  for  step  load  of  46,000  lbs  applied  at  10®  Ibs/sec 
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Figure  3.6.  MSHELL  model  w  displacement  vs.  time  for  step  load  of  46,000  lbs  applied  at  10® 
Ibs/sec 

In  the  next  case,  the  model  is  suddenly  loaded  at  a  rate  r  =  10*  Ibs/sec  in  an  effort  to  generate 
snap-through  in  both  directions  (snap-through  followed  by  snap-through  in  the  opposite  sense).  At 
this  higher  load  rate,  the  mass  does  snap  through  and  then  rebounds  completely  back  through 
the  w  =  d  plane.  Subsequent  oscillations  continue  to  cross  this  plane  (see  Figures  3.6  and  3.7V 
The  larger  amplitude  leads  to  a  lower  nonlinear  frequency  of  oscillation,  and  this  is  depicted  in 
Figure  3.8.  The  bifurcation  plot  for  this  load  rate  (not  shown)  looks  virtually  identical  to  that  of 
Figure  3.2,  except  that  the  onset  of  buckling  occurs  at  a  lower  load. 

S.S  Response  to  Sinusoidal  Loading 

To  simulate  the  behavior  of  a  node  in  the  DSHELL  finite  element  model,  in  which  eve  ry 
unrestrained  node  pulls  on  every  other  node,  the  MSHELL  model  is  loaded  sinusoidally.  The 
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Figure  3.7.  MSHELL  model  w  phase  diagram  for  step  load  of  46,000  lbs  applied  at  10*  Ibs/sec 


Frequency  (Hz) 

Figure  3.8.  MSHELL  model  frequency  response  for  step  load  of  46,000  lbs  applied  at  10*  Ibs/sec 
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Figure  3.9.  MSHELL  model  w  displacement  vs.  time  for  sinusoidal  load  of  10,000sin(1571t)  lbs 


loading  for  this  case  is 


1 


I  10,000sin(1571t) 


1 


(3.2) 


The  circular  frequency  of  1571  s“^  corresponds  to  /  =  250  Hz.  This  frequency  is  used  because 
it  approximates  the  nonlinear  fundamental  frequency  of  vibration  observed  in  the  finite-  element 
model  subjected  to  a  sudden  load  below  the  buckling  load  (see  Figure  4.13  on  page  4-12). 


Figure  3.9  plots  the  ui  displ^u;ement  of  the  mass  versus  time.  The  model  does  not  buckle 
until  the  fifth  cycle,  2rfter  which  snap  through  (in  both  directions)  occurs.  Though  not  shown  in 
Figure  3.9,  as  time  increases  the  model  is  seen  to  have  periods  of  oscillation  during  which  the  mass 
remains  on  the  prebuckled  side  of  the  w  =  d  plane  for  several  cycles,  then  resumes  the  buckled 
motion.  These  regions  are  seen  in  the  phase  diagram  (Figure  ;1.10)  as  small  orbits  near  the  value 
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Figure  3.10.  MSHELL  model  w  phase  diagram  for  sinusoidal  load  of  10,000sin(1571t)  lbs 

w  =  0  (Figure  3.10  represents  a  much  greater  time  span  than  Figure  3.9).  The  frequency  response 
of  this  model  (Figure  3.11)  shows  the  broad-b2U)ded  response  associated  with  chaos  (17:148). 

A  bifurcation  plot,  Figure  3.12,  was  generated  for  the  sinusoidal  loading  case.  In  the  region 
at  A,  the  orbits  suddenly  enlarge  —  the  result  of  a  flip  bifurcation.  The  origin  of  the  ‘hole’  at  B  is 
not  yet  understood,  but  the  points  C  represent  collapse  behavior. 

A  Poincare  map  of  the  case  with  amplitude  of  10,000  lbs  is  shown  in  Figure  3.13.  The 
Fourier  transform  indicates  chaotic  motion,  though  the  Poincare  map  shows  no  chaotic  attractor. 
This  is  due  to  the  lack  of  damping.  “If  a  system  does  not  have  sufficient  damping,  the  chaotic 
attractor  will  tend  to  fill  up  a  section  of  phase  space  uniformly  and  the  Cantor  set  structure,  which 
is  characteristic  of  strange  attractors,  will  not  be  evident  (17:134).” 
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Figure  3.11.  MSHELL  model  frequency  response  for  sinusoidal  load  of  10,000sin(1571<)  lbs  in¬ 
dicating  chaotic  behavior 
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Figure  3.12.  Bifurcation  plot  of  MSHELL  model,  sinusoidal  loading 
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Figure  3.13.  Poincare  map  for  sinusoidal  loading,  amplitude  of  10,000  lbs 


S.4  Methods  for  Choosing  Time  Step,  At 

Using  the  MSHELL  model,  a  series  of  simulations  was  run  to  find  etn  empirical  relationship 
between  the  system  dynamical  behavior  and  the  required  time  step.  At.  Two  cases  were  run  using 
the  MSHELL  model.  In  the  first  case,  the  MSHELL  model,  with  parameters  as  in  Table  3.1  on 
page  3-1,  is  loaded  sinusoidally  with  load  f  given  by 

1 
1 
1 

with  12  =  1571  radians/sec  (corresponding  to  /  =  250  Hz).  In  the  second  case,  the  model  is 
suddenly  loaded  at  a  load  rate  of  10®  Ibs/sec  to  a  load  fo  which  is  then  held  constant  for  all  time. 
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Figure  3.14.  MSHELL  model  (solid  line)  vs.  theory  (dotted  line)  for  =  2.6% 


The  load  function  is  therefore 


f  = 


1 


1  /(<). 


1 


where  f(t)  =  < 


10®  t 


fo, 


t  <  /o/lO® 
t  >  /o/lO® 


(3.4) 


For  each  value  of  /,  the  minimum  time  step  is  found.  The  criterion  for  a  successful  run  was  that 
the  maximum  cumulative  error,  Crmt,  never  exceeded  1%  over  the  course  of  the  run  of  0.2  seconds. 
The  1%  value  is  arbitrary,  but  Figures  3.14  and  3.15  indicate  the  improvement  gained  in  decreasing 
the  error  from  2.6%  in  Figure  3.14  to  1%  in  Figure  3.15. 

For  the  suddenly  applied  lozwl  of  infinite  duration,  there  is  a  clear  demarcation  between  the 
time  step  required  for  prebuckling  and  postbuckling  loads.  In  this  discussion,  postbuckling  implies 
a  reversal  of  geometry  about  the  supports  of  springs  3  and  4,  i.e.  the  mass  crosses  the  w  —  d 
plane.  Point  B  in  Figure  3.16  represents  the  lowest  load  at  which  buckling  occurred.  The  time 
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Figure  3.15.  MSHELL  model  (solid  line)  vs.  theory  (dotted  line)  for  Crm*  =  10% 

step  required  to  model  post  buckling  is  an  order  of  magnitude  less  than  that  required  to  model 
prebuckling.  Buckling  occurred  at  all  loads  beyond  that  of  point  B.  The  nonlinear  frequency  of 
oscillation  becomes  low  at  loads  just  prior  to  buckling,  and  this  is  the  reason  for  the  increase  in  At 
(point  A)  just  prior  to  the  buckling  load. 

For  the  sinusoidal  loading  case  there  is  no  clear  demarcation  between  pre-  and  postbuckling 
loads.  Over  the  course  of  the  0.2  second  run,  the  model  buckled  at  points  E  and  F  in  Figure  3.16, 
but  not  at  the  two  points  between  them.  Buckling  occurred  at  point  £  at  t  =  0.15  seconds,  and  at 
F  aX  t  =  0.042  seconds.  By  increasing  the  duration  of  the  run,  it  was  found  that  at  point  D  (the 
load  prior  to  point  E),  buckling  occurred  at  t  —  0.487  seconds.  Surprisingly,  at  a  still  lower  load 
(point  C),  buckling  occurred  at  (  =  0.236  seconds.  For  the  sinusoidal  loading  case,  the  buckling 
occurs  when  the  load  and  the  inertia  of  the  shell  model  are  close  enough  in  phase  for  a  long  enough 
time  to  allow  buckling.  These  conditions  occur  at  various  times,  and  a  higher  load  by  no  means 
assures  a  shorter  time  to  buckling. 
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Figure  3.16.  Maximum  time  step  At  vs.  load 


The  highest  natural  frequency  of  the  MSHELL  model  is  fh  =618  Hz,  corresponding  to  a  period 
of  T)i  =0.0016  seconds.  For  the  suddenly  loaded  case,  the  suitable  time  step  in  the  prebuckling  case 
is  about  t^/20,  while  for  the  postbuckling  case  a  time  step  of  about  t),/160  is  appropriate.  A  time 
step  in  this  region  lies  below  the  nearly  horizontal  line  containing  point  B  in  Figure  3.16  and  allows 
for  some  margin.  As  there  is  no  clear  demarcation  in  the  sinusoidal  loauiing  case,  the  suitable  time 
step  choice  is  still  about  71^/160. 
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IV.  THE  COMPOSITE  SHELL 


4.1  The  Composite 

A  24-ply  simply  supported  composite  cylindrical  is  the  subject  of  the  present  research.  The 
ply  lay-up  is  [Os/OOe],  and  the  material  properties  correspond  to  those  of  Hercules’  AS4-3501-6 
graphite  epoxy  composite.  The  material  properties  used  are  listed  in  Table  4.1.  Because  perfect 


wmsmi 

E2 

<J12 

Gi3 

C23 

U12 

18.844(10®) 

1.468(10®) 

0.91(10®) 

0.91(10®) 

0.45(10®) 

0.26 

Table  4.1.  Physical  properties  of  AS4-3501-6  graphite  epoxy  composite 


interply  bonding  is  assumed,  the  composite  is  modeled  as  a  4-ply  [0/90]j  lay-up  for  the  purposes 
of  input  to  the  finite  element  code.  The  geometry  and  loading  of  the  simply  supported  shell  are 
depicted  in  Figure  4.1.  The  first  three  natural  frequencies  are  found  from  linear  eigenvalue  analysis 
and  are  shown  in  Table  4.2.  The  edges  are  simply  supported  along  their  longitudinal  lengths,  and 


Table  4.2.  First  three  natural  frequencies  of  AS4-3501-6  graphite  epoxy  composite  shell 

the  sudden  load  is  applied  at  Node  1.  The  static  odlapse  load  of  the  shell  is  2800  lbs,  and  was 
calculated  using  the  method  of  Riks  (23).  The  shell  is  modeled  using  bilateral  symmetry,  which  is 
auiequate  for  cross-ply  laminates  (25:3-4).  Node  107  is  pointed  out  because  it  is  the  node  furthest 
from  DOF  restraints  due  to  either  boundaries  or  planes  of  symmetry  and  is  useful  in  the  analysis. 


4-2  Convergence 

Except  for  the  first  prebuckled  case  study,  the  mesh  used  in  the  current  research  is  the 
mesh  formulated  by  Taylor  (28:3-25).  Silva  (25:3-3)  performed  extensive  convergence  analyses  and 
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Figure  4.1.  Shell  geometry,  boundary  conditions,  and  loading  (Taylor’s  mesh) 

developed  criteria  for  guaranteeing  convergence  (see  Figure  4.2  for  Silva’s  mesh).  Taylor’s  mesh 
is  used  because  it  falls  within  five  percent  of  Silva’s  criteria  based  on  peak  load,  and  the  use  of 
Taylor’s  mesh  very  significantly  decreases  the  computer  run  time.  His  mesh  is  shown  in  Figure  4.3. 

^.S  Laminate  Failure 

In  the  present  research,  it  is  assumed  ply  failure  and  delamination  do  not  occur.  An  analysis 
by  Silva  (25:4-32)  using  the  two-dimensional  Tsai- Wu  failure  criterion  showed  that  ply  failure  begins 
very  early  with  the  initial  dimpling  of  the  shell,  and  become  extensive  when  rotations  exceed  2.5 
degrees.  He  notes,  however,  that  laminated  panels  have  been  observed  to  withstand  much  larger 
rotations  than  these  without  failure.  In  the  model  used  in  the  present  research,  energy  applied  to 
the  shell  in  the  form  of  the  external  load  can  go  into  elastic  strain  (potential)  energy  of  the  shell  or 
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plane  of  symmetry 


hinged  edge 


Figure  4.2.  Silva’s  mesh 
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Figure  4.3.  Taylor’s  mesh 
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free  edge 


shell  movement  (kinetic  energy)  alone.  Real  panels  have  additional  thermal  and  plastic  mechanisms 
to  dissipate  energy.  This  is  by  no  means  to  imply  that  the  snapping  shell  would  not  fail,  but  that 
classicEil  failure  criteria  applied  to  these  panel  predict  early  failure.  In  any  case,  in  order  to  study 
the  dynamic  snapping  phenomena  of  this  (and  other)  research,  the  ply  failure  is  ignored. 

4.4  Prebuckled  Behavior 

Prebuckled  cases  are  those  in  which  the  suddenly  applied  load  does  not  produce  collapse 
(buckling).  The  result  of  the  load  is  an  oscillation  of  the  shell  structure. 

4-41  The  Time  Step.  For  prebuckled  motion  of  the  shell  under  sudden  loading,  the  linear 
eigenvalue  analysis  may  be  used  to  formulate  the  appropriate  time  step.  The  time  step  may  be 
chosen  based  upon  the  first  linear  natural  frequency,  /i  : 

Ti  =  I//1 

At  =  ri/20  (4.1) 

This  time  step  was  found  to  be  sufficient  for  all  cases  of  prebuckled  motion. 

4-4-^  Case  Studies.  In  the  first  case,  a  convergence  tolerance,  etoi  of  0.2%  is  chosen  for 
the  equation  solver,  and  a  time  step  of  0.0001  seconds  is  used.  The  load  of  2000  pounds  is  appli' d 
at  a  rate  of  10®  Ibs/sec.  Silva’s  mesh  (25),  which  has  88  elements  (vs.  Taylor’s  49)  is  used.  The 
displacement  field  of  the  shell  structure  becomes  erratic  after  only  three  cycles  (see  Figure  4.1). 
Initially,  the  displacement  oscillates  about  the  node’s  static  solution  of  approximately  w  =  0.8.5 
inches,  but  then  the  displacements  become  large  and  erratic.  Furthermore,  these  large  erratic 
displacements  appear  in  coordinate  directions  other  than  that  of  the  applied  load  direction,  w.  In 
fact,  they  occur  in  every  degree  of  freedom,  as  can  be  seen  in  the  plots  of  displacements  at  tlie  free 
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Figure  4.4.  Displacement  of  Node  1,  the  point  of  applied  load,  e,<,i  =  0.2%,  Silva’s  mesh 

edge  node.  Node  163^  (see  Figure  4.5  and  Figure  4.6).  The  various  degrees  of  freedom  appear 
to  exhibit  instability  approximately  coincidentally  in  time  (at  a  particular  node).  Because  the  free 
circumferential  edge  of  the  shell  is  subject  to  the  least  constraint,  and  because  this  is  the  region 
in  which  the  in-plane  displacements  are  largest,  node  number  163^,  an  element-corner  node  on  the 
free  edge,  is  chosen  for  illustrative  purposes. 

Observing  the  motion  of  a  nodal  DOF  in  phase  space  (nodal  velocity  vs.  displacement)  is 
often  useful.  Consider  the  w  displacement  of  Node  1.  Figure  4.7  illustrates  the  nearly  elliptical 
trajectory  of  this  node.  This  is  the  center  node  of  the  model,  and  because  symmetry  has  been 
invoked,  all  DOF  at  this  node  have  been  constrained  except  it;.  The  very  regular  trajectory  is 
hence  a  by-product  of  shell  symmetry.  In  contrast,  consider  the  seven-DOF  node  163.  The  u,  i>, 
and  w  displacements  may  be  seen  in  Figure  4.5,  and  the  rotational  displacements  are  depicted  in 

'  Node  163  in  Silva’s  mesh  is  very  near  Node  107  in  Taylor'->  mesh. 

*Near  Node  107  in  Taylor’s  mesh 


4-6 


Figure  4.7.  Phase  space  of  w  displacement,  oscillating  region,  node  1,  e<o/  =  0.2%,  Silva's  mesh. 
Note  departure  of  the  trajectory  from  phase  orbit  due  to  numerical  instability 

Figure  4.6.  Tn  phase  space,  the  w  trajectory  in  the  oscillating  region  for  node  163  (see  Figure  4.8) 
is  much  less  regular  than  that  of  node  1.  The  trajectory  has  two  fairly  distinct  regions.  The 
roughly  circular  regions  of  the  trajectory  correspond  to  the  high  peaks  in  Figure  4.5,  while  the 
collection  of  points  near  the  origin  correspond  to  the  more  irreguleir  regions  between  those  peaks. 
The  unbounded  departure  of  the  trajectory  in  phase  space  indicates  numerical  instability,  as  the 
nodal  velocities  reach  clearly  unrealistic  values. 

Figures  4.9  and  4.10  depict  the  deformed  shell  geometries  prior  to  and  jdter  the  onset  of 
instability.  The  deformed  geometry  associated  with  instability  (Figure  4.10),  while  satisfying  com¬ 
patibility  and  potential  energy  constraints,  is  physically  unreasonable.  Some  adjacent  nodes  have 
collapsed  upon  each  other  to  form  odd  vertical  lines  and  many  elements  are  severely  missliapen. 
Examining  the  deformed  geometry  hence  provides  one  means  of  differentiating  between  dynamic 
and  numerical  instability. 
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Figure  4.8.  Phase  space  of  w  displacement,  oscillating  region,  node  163,  etoi  =  0.2%,  Silva’s  mesh 


P  =  2000  lbs 


0.002a  time 


Figure  4.9.  Deformed  shell  geometry,  numerically  stable  region,  etoi  =  0.2%,  Silva’s  mesh 


P  =  2000  lbs 


0.0028  time 


Figure  4.10.  Deformed  shell  geometry,  numerically  unstable  region,  etoi  =  0.2%,  Silva’s  mesh 

In  the  second  case,  two  changes  are  made.  First,  a  coarser  mesh  (Taylor’s)  is  used  and, 
second,  the  convergence  tolerance  for  the  Newton- Raphson  method  is  reduced  by  a  factor  of  2000 
to  0.0001%.  Both  of  these  changes  should  have  the  effect  of  improving  numerical  accuracy.  Reducing 
the  number  of  elements  reduces  the  number  of  nodes  by  nearly  a  half,  hence  reducing  the  number 
of  members  in  the  matrices  by  three-quarters.  The  time  step  is  unchanged.  The  results  indicate 
that  the  onset  of  instability  has  been  delayed,  as  no  instability  is  observed  over  the  length  of  the 
run. 

Figure  4.11  depicts  the  steady-state  oscillation  of  Node  1.  At  no  time  during  the  run  does  the 
displacement  diverge  as  it  did  in  the  previous  case.  The  use  of  the  tighter  tolerance  and  smaller 
time  step  hris  delayed  the  onset  on  instability  beyond  0.03  seconds.  Had  the  calculations  been 
continued  in  time,  an  instability  like  that  of  the  first  case  would  likely  have  been  seen.  In  fact, 
another  run  (not  included  here)  modelling  an  isotropic  shell  did  display  numerical  instability  just 
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Figure  4.11.  Displacement  of  Node  1,  e«o/  =  0.0001%,  Taylor’s  mesh 

one  time  step  prior  to  the  end  of  the  0.03  second  run.  The  phase  diagram  associated  with  activity  at 
Node  1  is  shown  in  Figure  4.12.  Figure  4.13  shows  the  fundamental  nonlineeir  frequency  to  be  241 
Hz.  This  nonlinear  frequency  is  roughly  half  that  of  the  fundamented  linear  frequency  of  vibration. 
Figure  4.14  shows  that  while  the  displacement  versus  time  at  the  shell  edge  is  periodic  eind  regular, 
the  function  appears  to  be  getting  noisier  (i.e.  more  high  frequency  content)  with  time.  It  is  not 
clear  whether  this  trend  is  due  to  parametric  resonance  or  the  onset  of  a  numerical  instability.  The 
trend  is  even  more  apparent  in  the  velocity  versus  time  (see  Figure  4.15). 

The  in-plane  motion  of  the  shell  is  substantial  even  in  this  pre-collapse  case  —  the  circum¬ 
ferential  (v)  displacements  are  about  one-third  of  the  transverse  {w)  displacements,  and  the  w 
displacements  are  about  eight  times  the  shell  thickness,  in  this  truly  nonlinear  problem,  the  as¬ 
sumption  of  in-plane  inextensibility,  common  in  linear  shell  analysis,  would  be  grossly  inadequate. 
The  phase  diagram  associated  with  activity  at  Node  107  is  shown  in  Figure  4.16. 
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Figure  4.13.  Frequency  response  of  Node  1,  etoi  =  0.0001%,  Taylor’s  mesh 
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Figure  4.14.  Displacement  of  Node  107,  eui  =  0.0001%,  Taylor’s  mesh 
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Figure  4.15.  Velocity  of  Node  107,  Ctoi  =  0.0001%,  Taylor’s  mesh 
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Figure  4.16.  Phase  trajectory  of  Node  107,  etoi  =  0.0001%,  Taylor’s  mesh 


The  frequency  response  at  the  shell  edge,  shown  in  Figure  4.17,  is  much  more  complex  than 
that  of  Node  1.  The  response  in  general  is  broad-banded  and  may  imply  a  near- chaotic  condition  or 
imminent  numerical  instability.  The  frequency  response  in  the  u  (longitudinal  in-plane)  direction 
of  the  entire  quarter  panel  of  the  shell  may  be  seen  in  Figures  4.18  and  4.19.  These  figures  represent 
the  relative  power  spectrum,  5(w),  for  the  fundamental  nonlinear  frequency  of  about  242  Hz  and 
the  third  harmonic  (near  763  Hz)  respectively.  In  the  figures,  S(w)  is  depicted  relative  to  the 
maximum  S{ui)  of  the  fundamental  frequency.  It  is  interesting  to  note  that  the  maocimum  values 
of  S{u)  in  Figure  4.19  do  not  lie  at  the  free  edge  of  the  shell,  but  rather  along  the  s  =  0  axis  of 
symmetry.  Similarly,  the  maximum  values  of  5(u;)  in  the  v  direction  take  place  along  the  x  =  0 
axis  of  symmetry. 

As  a  final  case,  an  isotropic  shell  of  the  same  physical  dimensions  was  analyzed.  The  elastic 
modulus  of  the  shell  was  adjusted  to  match  the  static  collapse  strength  of  the  composite  so  that 
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Figure  4.17.  Frequency  response  of  Node  107,  Ctoi  =  0.0001%,  Taylor’s  mesh 


the  effects  of  emisotropy  could  be  emphasized.  The  behavior  of  the  isotropic  shell  was  very  similar 
to  that  of  its  composite  counterpart.  However,  the  isotropic  shell  was  slightly  more  flexible  than 
the  composite.  The  maximum  displacements  were  about  5,  17,  and  6  percent  greater  in  the  u,  v, 
and  w  directions,  respectively,  for  the  isotropic  shell. 

4-5  Postbuckled  Behavior 

Postbuckled  behavior  refers  to  cases  involving  the  phenomenon  of  dynamic  buckling  or  snap- 
through  and,  in  this  context,  is  the  result  of  a  suddenly  applied  load  of  infinite  duration. 

The  methodology  for  choosing  the  appropriate  time  step(s)  is  an  integral  part  of  the  case 
studies  that  follow,  so  they  are  presented  together. 

In  Taylor’s  work  with  dynamic  cross-ply  shells,  the  smallest  time  step  used  was  At  =  0.0001 
seconds  (18:1652).  This  time  step  failed  to  capture  the  postbuckling  behavior  of  the  shell  (28:3-30). 
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Figure  4.18.  Power  Spectrum  5(w)  near  fundamental  nonlinear  frequency  of  oscillation,  /  =  242 
Hz,  for  u  displacements 
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Figure  4.19.  Power  Spectrum  S(w)  near  third  harmonic  of  fundamental  nonlinear  frequency  of 
oscillation,  /  =  763  Hz,  for  «  displacements 

As  a  first  attempt,  a  uniform  time  step  half  that  of  Taylor’s  (At  =  0.00005  seconds)  is  used  in 
evaluating  the  dyn2unic  snapping  of  the  cross-ply  shell.  The  w  displacement  vs.  time  is  seen  in 
Figure  4.20.  Prior  to  collapse,  the  highest  frequency  of  interest  is  not  associated  with  the  shell, 
but  with  the  suddenly  applied  load.  If  the  ramping  up  of  the  load  is  considered  to  be  the  first 
quarter  of  a  sinusoid,  the  ‘frequency’  of  the  applied  load  is  about  125  Hz,  The  time  step  of  0.00005 
seconds  corresponds  to  a  frequency  of  20,000  Hz  and  is  certainly  smaller  than  necessary  to  capture 
the  behavior  in  the  prebuckled  regime.  However,  the  time  step  does  not  appear  small  enough 
to  capture  the  behavior  after  buckling.  Figure  4.21  indicates  the  number  of  iterations  required 
to  converge  at  each  time  step.  The  Newton- Raphson  iteration  requires  an  inordinant  number  of 
iterations  to  converge  in  the  postbuckling  regime  (beginning  near  step  85),  at  one  point  requiring 
nearly  100  steps  to  converge.  An  examination  of  the  Node  1  velocity  data  indicates  that  the  run 
begins  to  produce  unrealistic  results  in  the  neighborhood  of  step  number  125,  as  the  nodal  velocities 
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Figure  4.20.  Displacement  vs.  time  for  fixed  At  =  0.00005  seconds,  Taylor’s  mesh 


Figure  4.21.  Number  of  iterations  required  vs.  time  step  for  fixed  At  =  0.00005  seconds 
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Figure  4.22.  Deformed  geometry  of  postbuckled  shell  indicating  numerical  instability,  At  = 
0.00005  seconds,  Taylor’s  mesh, 

become  extreme.  It  is  noteworthy  that  a  cursory  examination  of  the  nodal  displacement  versus  time 
(Figure  4.20)  does  not  give  an  indication  as  to  the  validity  of  the  data.  Once  again,  the  deformed 
geometry  (Figure  4.22)  gives  an  immediate  indication  that  something  has  gone  numerically  wrong. 
Figure  4.23  indicates  a  deformed  geometry  in  postbuckling  that  does  not  suffer  from  numerical 
instability. 

The  postbuckled  behavior  of  the  composite  cross-ply  shell  was  then  emalyzed  using  Taylor’s 
mesh  and  the  adaptive  time  step  method  described  in  section  3.4  on  page  3-11  with  otherwise 
identical  run  parameters.  A  time  step  based  upon  the  lowest  linear  frequency,  /j,  (from  the 
eigenvalue  analysis)  was  used  during  that  portion  of  the  run  describing  initial  collapse  of  the  shell. 
The  high-frequency  content  is  negligible  until  Node  1  reaches  its  meiximum  w  deflection,  i.e.,  the 
shell  ‘bottoms  out’.  Shortly  after  this  point  (in  this  run,  at  t  =  0.0049  sec)  parametric  (in-plane) 
resonances  cause  the  higher  frequencies  to  become  important.  Hence,  a  postbuckling  (collapsed) 
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0.002  s  time 


Figure  4.23.  Deformed  geometry  of  postbuckled  shell,  adaptive  time  step,  Taylor’s  mesh 

time  step  based  upon  the  highest  frequency  of  interest  seen  in  the  oscillating  shell,  /a  =  760  Hz, 
is  used  thereafter  (see  Figure  4.17  on  page  4-15),  The  highest  frequency  of  interest  was  arbitrarily- 
chosen  as  the  highest  frequency  within  10  dB  of  the  maximum  (from  data  at  the  most  flexible 
portion  of  the  structure,  i.e.  Node  107). 

Using  these  frequencies  with  the  criteria  developed  using  the  MSHELL  simplified  model,  we 

have 


"^pre 

=  l//i  ss  0.002sec 

^tpre 

=  Tpre/20  Si  O.OOOlsec 

(4.2) 

^po$t 

=  1//a  w  0.00131seconds 

Ripest 

=  rp„„/160  S5  8(10~®)seconds 

(4.3) 
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Figure  4.24.  Number  of  iterations  required  vs.  time  step  for  adaptive  At 


Table  4.3  indicates  the  result  of  using  the  adaptive  technique.  In  the  table,  Itot  is  the  total  number 


Method 

IKI 

CPU  seconds 

At  =  0.00005  (fixed) 

1698 

11.3 

20,000 

Adaptive  At 

2208 

5.9 

31,250 

At  =  Atpo,t  (est.) 

5500 

6 

80,000 

Table  4.3.  Comparison  of  time  step  schemes  for  buckling  of  composite  cross-ply  shell 


of  iterations  performed  over  the  run,  I  is  the  average  number  of  iterations  per  time  step,  and  CPU 
seconds  refers  to  CPU  seconds  of  run  time  on  the  CRAY  X-MP/216.  As  previously  mentioned,  the 
time  step  of  0.00005  seconds  failed  to  adequately  model  the  collapsing  shell  throughout  the  run. 
However,  Table  4.3  shows  that  using  Atpo,t  as  the  time  step  for  the  entire  run  is  not  a  good  option. 
The  adaptive  time  step  technique  offers  substantial  improvement.  Figure  4.24  indicates  how  the 
number  of  iterations  per  time  step  was  maintained  at  a  reasonable  level  throughout  tin-  run. 


As  in  other  schemes  for  nonlinear  analysis,  a  pnon  knowledge  of  the  sysi'  m's  ch.iracteristirs 
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Figure  4.25.  Displacement  vs.  time  of  Node  1.  adaptive  At,  Taylor’s  mesh,  collapse  case 


is  required.  However,  a  FFT  analysis  of  a  few  cycles  of  the  steady  state  oscillating  behavior  may 
provide  enough  data  to  proceed  with  minimum  trial  and  error.  The  time  step  for  this  trial  run  may 
be  chosen  based  on  the  linear  eigenvalue  analysis  as  in  Section  4.4  on  page  4-5. 

Figure  4.25  describes  the  w  displacement  Node  1,  and  Figure  4.26  is  the  phase  portrait  of 
the  motion  (for  the  postbuckled  region  only).  The  frequency  response  in  the  postbuckled  regime 
is  shown  in  Figure  4.27.  It  displays  quantitatively  the  same  sort  of  broad-banded  noise  associated 
with  the  chaotic  motion  of  the  MSHELL  simplified  model  (see  Figure  3.11  on  page  3-10). 

4-6  Propagation  of  Displacements 

The  nodes  in  Figure  4.28  were  chosen  to  attempt  to  observe  the  path  of  propagation  of 
displacements,  as  they  lie  approximately  along  a  line  from  Node  107  to  the  panel  center.  Figure  4.29 
indicates  how  the  amplitude  of  the  rebound  after  buckling  increases  moving  outward  from  the  panel 
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Figure  4.26.  Phase  diagram  of  Node  1,  Taylor’s  mesh,  collapsed  r'  gime  only 


Frequency  (Hz) 


Figure  4.27.  Frequency  response,  post  buckling,  adaptive  time  step 


Figure  4.28.  Shell  geometry  indicating  nodes  used  for  propagation  study 


center  and  is  greatest  at  the  free  edge.  The  three-dimensional  view  of  Figure  4.30  also  indicates  how 
the  displacement  at  the  edge  lags  the  center,  but  the  rebound  amplitude  is  larger.  The  frequency 
responses  of  these  nodes  is  shown  in  Figure  4.31.  All  nodes  display  the  broadband  frequency 
response  associated  with  chaos,  with  the  roll-off  (decrease  in  the  power  spectrum  with  increase  in 
frequency)  being  faster  further  from  the  panel  center. 
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Figure  4.31. 
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Frequency  response  of  nodes  lying  along  a  line  from  panel  center  to  panel  edge, 
postbuckling,  Taylor’s  mesh 
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V.  SUMMARY,  CONCLUSIONS,  AND  RECOMMENDATIONS 

5.1  Summary 

The  finite  element  model  used  in  this  research  (DSHELL)  incorporates  many  nonlinear  fea¬ 
tures  that  defy  closed  form  solution.  The  interaction  of  the  hundreds  i>l  displacement  and  rota¬ 
tional  degrees-of-freedom  complicate  analysis  of  the  problem.  At  the  present  time,  only  a  numerical 
(computer)  solution  is  available  for  solving  the  matrices  of  nonlinear  equations  involved.  Hence,  the 
numerical  parameters  take  on  major  importance.  Changes  in  the  time  step,  (At),  have  a  dramatic 
effect  on  the  problem  solution.  Time  steps  that  are  too  large  cause  the  unconditionally  stable 
(not  unconditionally  accurate)  integration  scheme  to  converge  to  an  incorrect  result  for  the  nodal 
displacements  and  rotations.  This,  in  turn,  results  in  erratic  and  unrealistic  deformed  geometries 
of  the  structure. 

In  an  attempt  to  study  the  features  of  this  complex  problem  on  a  smaller  scede,  a  simple  shell 
model  (MSHELL)  was  developed  capturing  many  of  the  features  of  the  DSHELL  finite  element 
model.  The  model  properties  were  adjusted  to  match  the  geometry,  static  collapse  load,  and  first 
linear  natural  frequency  of  the  finite  element  model.  Through  sinusoidal  excitation  of  the  model, 
the  mass  in  MSHELL  exhibits  behavior  in  many  respects  similar  to  that  of  a  single  node  in  the 
mesh  (excepting  rotational  behavior).  MSHELL  was  also  used  to  examine  the  sensitivity  of  the 
accuracy  in  the  numerical  integration  scheme  to  At. 

Much  of  the  nodal  displacement  data  were  taken  from  DSHELL  runs  and  transformed  through 
the  Fast  Fourier  Transform  (FFT)  to  the  frequency  domain.  The  FFT  d.tta  were  used  to  investigate 
amplitudes  of  the  vibrations  of  the  shell  at  various  frequencies.  .An  attempt  was  made  to  find 
trends  in  the  propagations  of  nodal  displacements  as  the  shell  vibrated,  and  a  comparison  was 
made  between  the  the  composite  shell  and  an  isotropic  shell  of  the  same  physical  dimensions  and 
static  collapse  strength. 
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5.2  Conclusions 


5.2.1  The  Time  Step  The  choice  of  time  step  for  the  finite  element  model  is  both  important 
and  difficult.  For  the  prebuckling  analysis,  computation  of  a  suitable  time  step  is  straightforward 
and  comes  from  the  linear  eigenvalue  analysis.  For  runs  involving  collapse,  an  adaptive  time  step 
method,  in  which  the  time  step  is  varied  over  the  course  of  the  run  in  accordance  with  the  criteria 
derived  from  the  MSHELL  simplified  model,  can  make  the  run  more  computationally  efficient. 
The  method  only  applies  to  cases  in  which  the  loading  frequency  (in  terms  of  load  rate,  for  the 
sudden  load  of  infinite  duration  case)  is  less  than  the  highest  frequency  of  interest  in  the  nonlinear 
oscillation  of  the  shell.  Where  this  is  not  true,  the  time  step  must  be  chosen  based  upon  the 
excitation  frequency. 

5.2.2  Chaos 

5.2.2. 1  Prebuckled  Case.  In  the  prebuckled  cases  examined,  chaoslike  behavior  was 
not  seen  either  in  the  MSHELL  model  or  the  DSHELL  finite  element  model.  Fourier  analyses  show 
that  despite  the  relatively  erratic  behavior  in  the  time  domain,  the  motion  is  a  combination  of 
oscillations  at  discrete  frequencies  associated  with  the  nonlinear  oscillation  of  the  system.  Frequency 
analysis  at  the  free  edge  node  did  indicate  a  general  broadening  of  the  response  indicating  near- 
chaotic  behavior  or  possibly  nearing  the  onset  of  numerical  instability.  In-p)ane  displacements  were 
found  to  be  significant  —  on  the  order  of  one-third  the  transverse  displacements.  This  parametric 
resonance  is  clearly  captured  by  DSHELL. 

Comparisons  between  the  composite  and  isotropic  shell  indicated  very  similar  behavior,  except 
that  the  maximum  displacements  in  the  isotropic  shell  (of  the  same  physical  dimensions  and  .static 
collapse  strength  as  the  composite  shell)  were  larger,  especially  in  ti,.-  circumferential  direction. 

5. 2. 2. 2  Postbuckled  Case.  Chaoslike  behavior  was  observed  in  both  the  MSHELL 
model  and  the  DSHELL  finite  element  model.  The  broad  banded  frequency  spectra  of  the  osril- 
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lations  and  the  very  irregular,  but  bounded,  trajectories  in  phase  space  are  consistent  with  chaos. 
Analysis  with  the  MSHELL  model  shows  that  the  dynamic  buckling  load  is  not  a  meaningful  pa¬ 
rameter  in  the  presence  of  cyclic  loading.  Buckling  occurs  under  specific  conditions  not  directly 
proportional  to  the  amplitude  of  the  forcing  function. 

5.2.3  Numerical  Instability  The  DSHELL  finite  element  code  has  a  threshold  for  numerical 
instability  that  is  iikely  a  function  of  the  following  parameters. 

•  the  number  of  unrestrained  DOF  (mesh  si*e) 

•  the  convergence  tolerance  of  the  Newton-Raphson  scheme 

•  the  total  number  of  iterations  performed  over  the  length  of  the  run  (a  function  of  the  time 
step  and,  to  a  lesser  extent,  the  convergence  tolerance) 

Cases  with  large  numbers  of  DOF  (fine  meshes),  large  convergence  tolerances,  and  large  time  steps 
will  suffer  numerical  instability  in  the  shortest  time.  Because  fine  convergence  tolerances  do  not 
significantly  increase  the  number  of  iterations  per  time  step,  a  small  tolerance  (e,,,/  <  10“®)  should 
routinely  be  used.  A  time  step  that  is  too  small  increases  the  number  of  calculations  unnecessarily 
and  hastens  the  onset  of  numerical  instability.  The  adaptive  time  step  method  presented  herein 
can  use  computer  time  more  efficiently. 

The  deformed  geometry  plots  give  an  immediate,  though  qualitative,  indication  of  numerical 
instability  while  observation  of  the  displacement  of  a  single  DOF  versus  time  may  not.  In  cases  of 
numericad  instability,  DSHELL  is  usually  able  to  converge  to  a  solution  that  satisfies  compatibility 
emd  potential  energy  constraints  while  arriving  at  an  unrealistic  deformed  geometry.  Unbounded 
trajectories  in  phase  space  also  indicate  numerical  instability  and  give  quantitative  information  ;i.> 
to  the  onset. 
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5.S  Recommendations 


5.3.1  The  Time  Step.  Unlike  multi-step  methods,  in  which  changing  the  time  step  At  is 
not  pragmatically  possible  during  the  run  (13:1346),  the  single-step  beta-m  method  can  easily 
accommodate  varying  At  (13:1350).  The  DSHELL  computer  code  could  be  modified  to  vary  the 
time  step  in  response  to  the  rate  of  change  of  one  or  more  variables.  In  this  way,  much  computer 
time  could  be  saved  by  lengthening  At  in  those  parts  of  the  run  (e  g.,  prebuckling  or  steady-state 
oscillation)  that  do  not  require  a  short  time  step.  This  would  also  allow  longer  runs  to  be  performed 
with  meaningful  results  (numerical  error  accumulation  could  be  much  delayed). 

5.3.2  The  Beta-m  Method.  The  beta-m  method  allows  the  user  to  modify  the  stability  and 
accuracy  of  the  algorithm  through  the  choice  of  the  coefficients,  fii.  For  the  current  (and  prior) 
research,  the  /?,-  were  chosen  to  yield  unconditional  stability.  Optimum  choices  of  these  coefficients 
may  exist  based  upon  the  type  of  problem  being  investigated.  Confidence  gained  in  the  choice 
of  time  step  allows  the  investigator  to  shift  emphasis  from  unconditional  convergence  to  higher 
accuracy. 

5.3.3  Convergence  Tolerance.  In  the  current  research,  the  convergence  tolerance  used  for 
the  Newton-Raphson  iteration  technique  weis  practically  eliminated  as  a  variable  by  choosing  a 
minute  value.  Relaxation  of  the  convergence  tolerance  might  yield  somewhat  faster  computer  runs, 
though  the  interaction  between  the  convergence  tolerance  and  the  time  step  is  not  fully  understood. 
The  time  step  has  a  more  pronounced  effect  on  the  number  of  iterations  required  for  convergence. 
Future  reseuch  could  investigate  this  relationship. 

5.3.4  Modem  Nonlinear  Analysis.  Modern  nonlinear  analysis  involves  a  host  of  techniques 
not  applied  in  the  current  research.  These  include  study  of  the  Lyapunov  exponents,  fractal  dimen¬ 
sions,  and  many  mor  graphical  techniques  (logistic  maps,  circle  maps,  horseshoe  maps,  etc.)  not 
presented  here.  While  the  DSHELL  model  would  be  too  cumbersome  for  most  of  these  methods, 
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the  MSHELL  model  could  be  used  to  explore  many  of  the  phenomena  associated  with  buckling 
shells. 

5.3.5  Optimization  of  DSHELL.  The  DSHELL  code  could  be  optimized  to  choose  an  initial 
time  step  and  convergence  tolerance  based  on  the  problem  presentation.  The  program  could  adjust 
these  parameters  based  upon  the  dynamic  behavior  involved.  For  example,  the  program  could  be 
made  to  recognize  the  counter-flexural  point  associated  with  imminent  snap-buckling  and  adjust 
the  time  step  accordingly.  Vectorization  of  the  current  DSHELL  code  is  probably  not  practical, 
but  a  new  vectorized  code  based  upon  the  same  theory  could  be  written.  This  would  be  a  major 
tcisk,  but  the  time  involved  in  running  the  current  code  (even  on  the  supercomputer)  is  formidable. 


Appendix  A.  Derivation  of  MSHELL  Equations  of  Motion 

The  simplified  shell  model  is  depicted  in  Figure  2.3  on  page  2-9,  and  its  appearance  as  it 

relates  to  the  shell  geometry  is  shown  in  Figure  2.4  on  page  2-10.  All  springs  are  of  undeformed 

length  6,  and  springs  1  and  2  remain  horizontal  regardless  of  the  deformed  geometry  (the  ends  of 

springs  1  and  2  can  be  thought  of  as  riding  in  tracks  that  have  no  resistance  to  vertical  motion 

of  the  spring  ends).  The  object  is  of  mass  m,  and  the  dimensions  h  and  d  are  the  projections 

of  undeformed  spring  3  (or  4)  on  the  horizontal  and  vertical,  respectively.  Given  the  vector  of 

T 

displacements  x  =  [  y  y  u;  ]  •  deformed  spring  lengths,  are 

t>i  =  [(b  -  uf  + 

63  =  [(6-hu)2-H,;2]^/" 

63  =  [(/» -t- +  (d  -  tt;)^ -h  u^] 

64  =  [(/i  -  -I- (d- U))^ -I- (A.l) 

and  the  change  in  length  of  each  spring  is  given  by  A;  =  6,  —  6; 

Ai  =  5j  _  6  =  [(6  -  u)2  4- -  6 

Aj  =  b2-b=  [(b-f-uf  +  v^]''^ -6 

A3  =  63  -  6  =  [(ft  4- v)- -I- (d  -  u>)^ -I- u'] -  6 

A4  =  64  —  6  =  [(A  —  4- (d  —  u>)^  4- u^]  —  6  (-^.2) 

The  force  generated  by  each  of  the  springs  is 

fi=l-,Aie,  (A.;)) 


A- 1 


where  ki  are  the  linear  spring  stiffnesses.  The  direction  of  each  of  liie  forces  is  given  by 


V. 

IIV.II 


(A.4) 


where  V,  represents  a  vector  of  length  bi  pointing  along  the  spring  in  the  direction  of  the  force 
generated  by  a  positive  Ai  .  The  unit  vectors,  Ci  ,  are  thus  easily  written  component-wise  as 


ei  =  [(6- u)i- iij-)-0k]/6i 

eo  =  [—(6  -(-  u)i  —  ij  -t-  OkJ/feo 
©3  =  [-«i  -  (h  +  elj  +  (d  -  ie)k]/63 

64  =  [-ui  +  (h  -  v)j  +  {d- w)k]/b4  (A. 5) 

Summing  the  forces  contributed  by  the  deformed  springs,  we  have 


f  =  ^fciAie.  (A.6) 

I 

We  include  a  very  simple  model  of  viscous  damping  by  assuming  that  the  velocity  of  the  mass  i.s 
resisted  equally  in  all  directions  by  an  amount  proportional  to  the  velocity  in  that  direction.  That 
is 

^damj'in)  —  ~CX  (.A.i  ) 

where  c  >  0  is  the  damping  coefficient  and  x  represents  the  vector  of  velocities  ii,  ii,  and  iv.  This  is 
a  significant  simplification  over  the  DSHELL  model,  which  utilizes  a  proportional  damping  scheme: 


C  =  2^-M 


(A.8) 


.\-2 


Where  ^  is  the  damping  ratio  and  w  is  the  linear  fundamental  natural  frequency  (6:727).  We  further 
allow  for  a  forcing  function  that  may  vary  with  time: 


f /arcing  —  f(t) 

Setting  these  forces  equal  to  the  inertial  force  of  the  mass  yields 

“  ^springs  "h  ^damping  "h  ^forcing 

—  ki^ie[  T  k2^2^2  +  k^^^e^  +  1:4^404  —  cx  +  f(t) 

Substituting  Eqs  (A. 2)  and  (A. 5)  into  Eq  (A.  10)  yields 


mx  = 


fjprinjj  "h  ^damping  "t"  ^forcing 

=  ki{bi  -  b)  ^(6  -  u)i  -  tj  +  Ok  /6i  + 

+k2{b2  -  b)  |-(6  +  u)i  -  uj  +  Ok  /b^  + 
+*;3(*3  -  b)  |-ui  -  (h  +  t;)j  +  (d  -  u  )kj  763  + 
+^4(64  -  b)  J-«i  +(h-  i;)j  +  (d  -  u;)k|  /b^  - 
— cx  +  f(t) 


which  may  be  factored  as 


mx  =  1 


*1  (6  -  u)  - 
bk^  (6  +  u) 


bk\  {b  —  u) 


62 


+  j 


bi 

-CU  +  fu 

,  6*4 {h  - 

-  v) 

64 

6*3  (h  +  f ) 

bksu  bkgu  ,  ,, 

— - kgU-\ - 7 - ^2  (^  +  u) 

t>3  64 


,  bk\v  ,  bkav 
-kiv+  -^-k2V+  — ^ 
0i  62 


bs 


A-3 


(A.9) 


(A.IO; 


(A.ll) 


,  .  bk3(d-w)  bk^id-w) 

ks  {d-  w) - 1 - -  +  k^id  -  w)  - 


+  fc 

-CW  +  fui] 


63 


64 


(A. 12) 


Gathering  by  vector  component  yields 


»  bki(b~u)  bksu  bk^u 

1:  mu  =  fci(6-u) - 7 - K3U+— - K4U-I — - — 

0%  f>3  Oa 

-  to  (6  +  u)  +  -CU  +  fu 


j  ;  mv  =  k4(h  —  v)  — 


f>2 

bk4  {h  —  y) 


,  bk\v  ,  bk-iV 

-kiv+— - k2V+-— 

O4  Oi  02 


k  :  mii;  = 


-  k3{h  +  v)+  _  cy  ^ 


64 


-  CW  +  fvj 


(A.13) 


These  equations  yield  closed  form  solutions  for  the  three  acceleration  components  in  terms  of 
the  displacement,  velocity,  and  forcing  function  components.  These  equations  (factored  slightly 
differently)  are 


u(u,y,u),u,/„(t))  = 

(1/m)  +6^1 - -  -  bk2  +  —~ 

v(u,v,w,v,fv{t))  = 

/,  /  \  . .  bhkz  ,  ,  bhk^ 

(1/m)  j/„  -  hkz  +  +  /i*4  - 

+  [-ki  +  -ir-k2+-rr-k3  +  -^-k4+-r^]v-t 
\  61  62  63  64  J 

ti)(u,  V,  w,  w,  fu,(t))  = 


(A.14) 


(1/m) 


/u,  +  dk^ - 7 - 1-  dk^  —  —7 —  (  “^3  ”7 - M  H — ; —  )  ^  “ 

03  64  V  63  64  / 


cw 


(A. 15) 


(A. 16) 
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These  equations  may  be  written  in  the  form 


[M]{x}  +  [C]{x}  +  {[K],,„,<.,  +  [K]  nonlinear'^  {^1  {f}/orctnj 


Which  in  expanded  form  is 


where 
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Pi 
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^1  +  ^2  +  ^'3  +  ^4 
^3  +  ^4 


=  6 


6  (—— 

\  6i  62  63  64  u  u  y  V^i“  ^2“ 

/  ^1  ^2  ^3  ^4  A^3  hk4\  hkz  hk^ 

\  61  62  ^3  ^4  bsV  b^v )  V  V 


^2  ^3  ^4  A^3 

62  63  64 

<>*4  ^  J  ^  ^  ^  d<r3  dt4 

64  \  63  63U)  64U;/  w  w 


amd  the  6,  are  found  in  Eq  (A.l). 
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(A.l?) 


(A.18) 


Appendix  B.  PREDCORR:  Prcdictor-Corvtctor  Subroutine 


This  FORTRAN  subroutine  implements  the  predictor-corrector  algorithm  described  in  Section  2.2. 
on  page  2-15. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  predcorr 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
iaplicit  double  precision  (a-h,o-z) 

conion  /uT«/  uO,TO,BO,udotO,TdotO,sdotO,u,v,*,udot,Tdot,edot, 
k  uddot , vddot , vddot 

couon  /phys/  spm^kl,sprngk2,spmgk3,8prngk4,iiiass,h,b,d,c_daap 
couon  /forces/  f_a*pl,f_freq,f_dura,fu_dir,fT_dir,fB_dir ,fu,f», 
ft  fw,force_u,force_»,force_B,nftype 
conaon  /tiae/  del_t,t_stop,tia, epsilon, r_load,t_raap,del_t_last, 
ft  del_t _aax , del_t 1 , chg_t 1 , del_t2 , chg_t2 ,snap_t iae , t_shif t , if ix_t 

coaaon  /iterate/  r_iter _tot, iter, j_8tep, iter _aax,i_snap_f lag 
coaaon  /aisc/  pi,rad8,dega 

coaaon  /printer/  i_print,idat_print,ipm_intrTl,idat_intrvl, 
ft  i_3nap ,  j_print ,  jdat_print ,  i.screen,  i_auto_try ,  ipmO ,  ipdatO ,  if  ft 
coaaon  /old/  u_laat,v_la8t  ,B.laat,fla8t_u,fla8t_»,flast_B,B_tot, 
ft  udot.last ,Tdot_la8t,Bdot_last,uddot_last ,Tddot_last ,Bddot_last 
coaaon  /error/  err_total,pct_err_a2ut,pct_error,top_err , 
ft  top_err_tia,erT_cua,ra8_err,pct_err_8ua,i_bif_Brite,i.tcheck 
coaaon  /energy/  bl,b2,b3,b4,pe_l,pe_2,pe_3,pe_4,e_potential, 
ft  e_kinetic,e_pot_la8t,e_kin_last 
coaaon  /beta/  beta_0,beta_l,beta_2,aethod 

coaaon  /univ/  i_on_unT,i.unT_iiit,i_disp,i_Tel,i_acc,icase,ncount 
coaaon  /poinc/  p_pha8e,tau,i_on_poin,i_poiii_u,i_poin_v,i_poin_B 
coaaon  /bifur/  force_8top,i_on_bifur,n_force8,n_ignore,n_capture, 
ft  i_bif _8tep , i_f orces , i_bif _u , i_bif _v , i.bif _b , i.err.f lag , i_good 

c 

c  j.step  •  tiae  step 

c  iter  ■  iteration  step  eithin  tiae  step  j_step 
c 

iter  “  1 

c  write (8,*) ’iter  “  ’,iter 

c 

c  for  first  iteration,  estiaate  velocity  at  end  of  jtb  tiae  step, 
c 

udot.guess  »  udot  >  uddot*del_t 

u_gue8s  "  u  +  0.6*udot*del_t  +  0 . 5*udot_gue88*del_t 
c 

vdot.guess  >  Tdot  >  Tddot«del_t 

▼-guess  ■  B  +  O.Seydotedel.t  +  0 . 6*Tdot_gues8*del_t 
c 

Bdot.guess  >  Bdot  -f  Bddot«del_t 

B.guess  «  B  +  0.5*Bdot*del_t  +  0 . 5*Bdot_gue88*del_t 
c 

call  a.calc (u-guess , v.guess , B-guess , udot .guess ,vdot .guess , 
ft  Bdot.guess ,uddot.guess,vddot.guess,Bddot.guess) 

c 

c  for  iterations  greater  th^ul  first _ 

c 

20  udot.guess  ■  udot  +  0.5*uddot*del_t  +  0.5*uddot.gues8*del.t 


B-1 


ulast  ■  u_guass 

u.guess  ■  u  +  0.5*udot*del_t  +  0 . 5*udot_gue88*del-_t 
c 

vdot.guess  ■  Tdot  +  0 . 5*Tddot*del_t  +  0 . 5*vddot_gue8s*del_t 
▼last  ■  T_gues8 

v_gue88  »  V ■+  0.5*Tdot*del_t  +  0.5*7dot_guea8*del_t 
c 

vdot_gue88  =  Bdot  +  O.S*wddot*del_t  +  0.5*wddot_gue8s*del_t 
Blast  »  B.guess 

B_gues8  «  B  +  0.5*Bdot*del_t  +  0.5*Bdot_gu«88*del_t 
c 

call  a_calc (u.guess , v.gueas , B.gueas ,udot_guesa , vdot .guess , 
k  Bdot .guess, uddot .guess, Tddot .guess, Bddot.guess) 

c 

c  check  for  convergence 
c 

30  ep8chk2  =  dsqrtC  u.guesa**2  +  v.guea8**2  +  B.guess**2) 
epschkl  «  dsqrt(  ulast**2  +  vla8t**2  +  Bla8t**2) 
epschk3  •  dab8(epschk2  -  epschkl) 
epschk4  >  ep8ilon*epschk2 
if  (ep8Chk3  .gt.  epachk4)  then 

if  (i.print  .eq.  1  .and.  j.print  .eq.  iprn.intrvl) 
k  Brite(8,*) ’iter  ••  ’,iter,’  epschk3  =  ’,epschk3 

if  (iter  .gt.  iter .max)  then 
if  (i.auto.try  .ne.  1)  then 
call  report.it 
stop 
end  if 

i.err.flag  ■  2 
call  err.hndlr 
return 
end  if 

c  BTite (20 , 1000) j.step, iter , B_gues8 

iter  •  iter  +  1 
r.iter.tot  "  r.iter.tot  +  1.0 
goto  20 
end  if 

40  if  (i.print  .eq.  1  .and.  j.print  .eq.  iprn.intrvl) 
k  BTite (8, ♦) ’iter  “  ’,iter,’  epschk3  •  ’,ep8chk3 

u. last  ■  u 

v. last  “  V 
B.last  >  B 
u  *  u.guess 
▼  ■  v.guess 
*  ~  B.guess 

udot  ■  udot.guesB 
vdot  >  vdot .guess 
Bdot  -  Bdot .guess 
uddot  s  uddot.guess 
vddot  >  vddot.guess 
Bddot  *  Bddot.guess 
1000  forBat(lx,i5,lx,i4,lz,el2.5) 
return 
end 


B-2 


Appendix  C.  CHECK.IT:  Energy  Checking  Subroutine 

This  FORTRAN  subroutine  implements  the  energy  checking  scheme  described  in  Section  2.2.4  on 
page  2-13. 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
subroutine  check. it 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
implicit  double  precision  (a-h,o-z) 

common  /uve/  uO,vO,eO,udotO,vdotO,wdotO,u,v,s,udot,vdot,Bdot, 
ft  uddot , Tddot , Bddot 

common  /phys/  spmgkl,sprngk2,spmgk3,spmgk4,rmas8,h,b,d,c_damp 
common  /forces/  f_ampl,f_freq,l_dura,fu_dir,f»_dir,fH_dir ,fu,fv, 
ft  lB,force_u,force_T,force_H,nftype 
common  /time/  del_t,t_stop,tim, epsilon, r_load,t_ramp,del_t_last, 
ft  del_t_max,del_tl ,chg_tl ,del_t2 ,chg_t2,snap_time, t .shift , if ix_t 

common  /iterate/  r.iter.tot , iter, j. step, iter .max, i.snap.f lag 
common  /misc/  pi,rads,degs 

common  /printer/  i.print,idat_print,ipm_intrvl,idat.intrvl, 
ft  i.snap , j. print , jdat.print , i.screen, i.auto.try , iprnO , ipdatO , if f t 

common  /old/  u.la8t,v.last,B.la8t,fla8t.u,fla8t.v,flast_B,B.tot, 
ft  udot. last ,vdot. last (Bdot.last ,uddot.last ,Tddot.last ,Bddot.last 

common  /error/  err.total,pct.err.max,pct.error,top.err , 
ft  top.err.t im,  err.cum ,rm8.err , pct.err.sum , i.bif .bt ite , i.tcheck 

common  /energy/  bl,b2,b3,b4,pe.l,pe.2,pe.3,pe.4,e.potential, 
ft  e.kinetic,e_pot_la8t,e.kin.last 
common  /beta/  beta_0,beta_l,beta_2, method 

common  /univ/  i_on_unT,i.unv.int,i_di8p,i.Bel,i.acc,icase,ncount 
common  /poinc/  p.pha8e,tau,i.on.poin,i.poin_u,i_poin.T,i_poin.B 
common  /bifur/  force_8tep,i.on.bifur,n.forceB,n.ignore,n.capture, 

*  i.bif .step , i.f orces , i.bif  _u , i.bif .▼ , i.bif .b , i.err.f lag , i.good 

c 

c  calculate  potential  energy  of  system  by  spring 
c 

call  b.calc(u,T,B) 
c 

pe.l  =  0.5*8pmgkl*(  bl  -  b  )**2 
pe_2  “  0.5*8pmgk2«(  b2  -  b  )**2 
pe_3  ■  0.5espmgk3*(  b3  -  b  )**2 
pe.4  ■  0.5*8pmgk4*(  b4  -  b  )**2 
c 

s.potantial  •  pe_l  +  pe.2  +  pe.3  +  pe.4 
c 

c  calculate  kinetic  energy  of  mass 

c 

e.kinetic  ■  O.Sermass  *  (  udot**2  Bdot**2  Bdot**2  ) 
c 

c  calculate  change  in  energies  over  tine  step 
c 

del.e.pot  ■  e.potential  -  e.pot.last 
del.e.kin  ■  e.kinetic  -  e.kin.last 
delta.e  ■  del.e.pot  +  del.e.kin 
c 

c  calculate  increment  of  Bork  done  by  applied  load 
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c  calculate  displaceaent  over  last  increaent 
c 

dal_u  ”  u  -  u.laat 
del_v  ■  T  -  T.last 
del_B  ■  w  -  B.laat 

c  writedO,*) ’del_? — >u«’,del_u,’  v«’,del_v,’  B*’,del_B 

c 

c  calculate  Bork  done  over  last  increment  using  avg  force  over  step 
c 

Bork_u  =  0.5*(force_u  +  flast.u)  •  del_u 
Bork_v  «  0.5*(force_v  +  flast_v)  •  del_v 
Bork_B  •  0.5*(force_B  +  flast_B)  •  del_B 
del_Bork  »  Bork_u  +  Bork_v  +  Bork_B 
c  BriteClO,*) ’Bork_? — >  u*‘ ,Bork_u, ’  v«’,Bork_v,’  B=’,Bork_v 

c 

c  add  to  total  Bork  done 
c 

B_tot  =  B_tot  +  del_Bork 
c 

c  take  sum  (should  remain  constant  thoughout  run  in  undamped  case) 
c 

err_3tep  »  daba(del_Bork  -  del_e_pot  -  del_e_kin) 
err_cum  »  err_cum  +  err.step 
c 

c  calculate  percent  error 
c 

if(B_tot  .gt.  l,0D-10)pct .error  ■  dabs (err.cum/B.tot)* 100. 
pct.err.sum  ■  pct.err.sum  +  pct_error**2 
rms.err  =  dsqrtCpct.err.sum/j.step) 
if  (rms.err  .gt.  top.err)  then 
top.err  »  rms.err 
top.err.tim  ■  tim 
end  if 
c 

c  check  for  exceeding  pct.err.mcix 
c 

if  (rms.err  .gt.  pct.err.max)  then 
i.err.flag  -  1 
if  (i.auto.try  .ne.  Dthen 
call  report.it 
stop 
end  if 

call  err.hndlr 
return 
end  if 
c 

c  BritedO,*) ’B.tot"  e.kinetic"  ’.e.kinetic 

c  BritedO,*) ’e.potential*  ’  .e.potential,  ’  err_total»  ’.err.total 

return 
end 
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